laura is hosted by Hepforge, IPPP Durham
Laura++  v3r2
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauBelleCPCoeffSet.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2006 - 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 #include <fstream>
17 #include <vector>
18 
19 #include "TMath.h"
20 #include "TRandom.h"
21 
22 #include "LauComplex.hh"
23 #include "LauConstants.hh"
24 #include "LauBelleCPCoeffSet.hh"
25 #include "LauParameter.hh"
26 #include "LauPrint.hh"
27 
29 
30 
31 LauBelleCPCoeffSet::LauBelleCPCoeffSet(const TString& compName, Double_t a, Double_t delta, Double_t b, Double_t phi,
32  Bool_t aFixed, Bool_t deltaFixed, Bool_t bFixed, Bool_t phiFixed, Bool_t bSecondStage, Bool_t phiSecondStage) :
33  LauAbsCoeffSet(compName),
34  a_(new LauParameter("A", a, minMagnitude_, maxMagnitude_, aFixed)),
35  b_(new LauParameter("B", b, minMagnitude_, maxMagnitude_, bFixed)),
36  delta_(new LauParameter("Delta", delta, minPhase_, maxPhase_, deltaFixed)),
37  phi_(new LauParameter("Phi", phi, minPhase_, maxPhase_, phiFixed)),
38  particleCoeff_(0.0,0.0),
39  antiparticleCoeff_(0.0,0.0),
40  acp_("ACP", (-2.0*b*TMath::Cos(phi))/(1.0+b*b), -1.0, 1.0, bFixed&&phiFixed)
41 {
42  if (bSecondStage && !bFixed) {
43  b_->secondStage(kTRUE);
44  b_->initValue(0.0);
45  }
46  if (phiSecondStage && !phiFixed) {
47  phi_->secondStage(kTRUE);
48  phi_->initValue(0.0);
49  }
50 }
51 
52 LauBelleCPCoeffSet::LauBelleCPCoeffSet(const LauBelleCPCoeffSet& rhs, CloneOption cloneOption, Double_t constFactor) : LauAbsCoeffSet(rhs.name()),
53  a_(0),
54  b_(0),
55  delta_(0),
56  phi_(0),
57  particleCoeff_( rhs.particleCoeff_ ),
58  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
59  acp_( rhs.acp_ )
60 {
61  if ( cloneOption == All || cloneOption == TieMagnitude ) {
62  a_ = rhs.a_->createClone(constFactor);
63  } else {
64  a_ = new LauParameter("A", rhs.a_->value(), minMagnitude_, maxMagnitude_, rhs.a_->fixed());
65  if ( rhs.a_->blind() ) {
66  const LauBlind* blinder = rhs.a_->blinder();
67  a_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
68  }
69  }
70 
71  if ( cloneOption == All || cloneOption == TieCPPars ) {
72  b_ = rhs.b_->createClone(constFactor);
73  } else {
74  b_ = new LauParameter("B", rhs.b_->value(), minMagnitude_, maxMagnitude_, rhs.b_->fixed());
75  if ( rhs.b_->blind() ) {
76  const LauBlind* blinder = rhs.b_->blinder();
77  b_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
78  }
79  }
80 
81  if ( cloneOption == All || cloneOption == TiePhase ) {
82  delta_ = rhs.delta_->createClone(constFactor);
83  } else {
84  delta_ = new LauParameter("Delta", rhs.delta_->value(), minPhase_, maxPhase_, rhs.delta_->fixed());
85  if ( rhs.delta_->blind() ) {
86  const LauBlind* blinder = rhs.delta_->blinder();
87  delta_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
88  }
89  }
90 
91  if ( cloneOption == All || cloneOption == TieCPPars ) {
92  phi_ = rhs.phi_->createClone(constFactor);
93  } else {
94  phi_ = new LauParameter("Phi", rhs.phi_->value(), minPhase_, maxPhase_, rhs.phi_->fixed());
95  if ( rhs.phi_->blind() ) {
96  const LauBlind* blinder = rhs.phi_->blinder();
97  phi_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
98  }
99  }
100 }
101 
102 std::vector<LauParameter*> LauBelleCPCoeffSet::getParameters()
103 {
104  std::vector<LauParameter*> pars;
105  pars.push_back(a_);
106  pars.push_back(b_);
107  pars.push_back(delta_);
108  pars.push_back(phi_);
109  return pars;
110 }
111 
113 {
114  std::cout<<"INFO in LauBelleCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
115  std::cout<<"a-magnitude = "<<a_->value()<<",\t";
116  std::cout<<"delta = "<<delta_->value()<<",\t";
117  std::cout<<"b-magnitude = "<<b_->value()<<",\t";
118  std::cout<<"phi = "<<phi_->value()<<"."<<std::endl;
119 }
120 
121 void LauBelleCPCoeffSet::printTableHeading(std::ostream& stream) const
122 {
123  stream<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
124  stream<<"\\hline"<<std::endl;
125  stream<<"Component & a-Magnitude & delta & b-Magnitude & phi \\\\"<<std::endl;
126  stream<<"\\hline"<<std::endl;
127 }
128 
129 void LauBelleCPCoeffSet::printTableRow(std::ostream& stream) const
130 {
131  LauPrint print;
132  TString resName = this->name();
133  resName = resName.ReplaceAll("_", "\\_");
134  stream<<resName<<" & $";
135  print.printFormat(stream, a_->value());
136  stream<<" \\pm ";
137  print.printFormat(stream, a_->error());
138  stream<<"$ & $";
139  print.printFormat(stream, delta_->value());
140  stream<<" \\pm ";
141  print.printFormat(stream, delta_->error());
142  stream<<"$ & $";
143  print.printFormat(stream, b_->value());
144  stream<<" \\pm ";
145  print.printFormat(stream, b_->error());
146  stream<<"$ & $";
147  print.printFormat(stream, phi_->value());
148  stream<<" \\pm ";
149  print.printFormat(stream, phi_->error());
150  stream<<"$ \\\\"<<std::endl;
151 }
152 
154 {
155  if (a_->fixed() == kFALSE) {
156  // Choose an a-magnitude between 0.0 and 2.0
157  Double_t mag = LauAbsCoeffSet::getRandomiser()->Rndm()*2.0;
158  a_->initValue(mag); a_->value(mag);
159  }
160  if (b_->fixed() == kFALSE && b_->secondStage() == kFALSE) {
161  // Choose a b-magnitude between 0.0 and 0.1
162  Double_t mag = LauAbsCoeffSet::getRandomiser()->Rndm()*0.1;
163  b_->initValue(mag); b_->value(mag);
164  }
165  if (delta_->fixed() == kFALSE) {
166  // Choose a phase between +- pi
168  delta_->initValue(phase); delta_->value(phase);
169  }
170  if (phi_->fixed() == kFALSE && phi_->secondStage() == kFALSE) {
171  // Choose a phase between +- pi
173  phi_->initValue(phase); phi_->value(phase);
174  }
175 }
176 
178 {
179  // retrieve the current values from the parameters
180  Double_t aVal = a_->value();
181  Double_t bVal = b_->value();
182  Double_t deltaVal = delta_->value();
183  Double_t phiVal = phi_->value();
184  Double_t genDelta = delta_->genValue();
185  Double_t genPhi = phi_->genValue();
186 
187  // Check whether we have a negative "a" magnitude.
188  // If so make it positive and add pi to the "delta" phase.
189  if (aVal < 0.0) {
190  aVal *= -1.0;
191  deltaVal += LauConstants::pi;
192  }
193 
194  // Check whether we have a negative "b" magnitude.
195  // If so make it positive and add pi to the "phi" phase.
196  if (bVal < 0.0) {
197  bVal *= -1.0;
198  phiVal += LauConstants::pi;
199  }
200 
201  // Check now whether the phases lies in the right range (-pi to pi).
202  Bool_t deltaWithinRange(kFALSE);
203  Bool_t phiWithinRange(kFALSE);
204  while (deltaWithinRange == kFALSE && phiWithinRange == kFALSE) {
205  if (deltaVal > -LauConstants::pi && deltaVal < LauConstants::pi) {
206  deltaWithinRange = kTRUE;
207  } else {
208  // Not within the specified range
209  if (deltaVal > LauConstants::pi) {
210  deltaVal -= LauConstants::twoPi;
211  } else if (deltaVal < -LauConstants::pi) {
212  deltaVal += LauConstants::twoPi;
213  }
214  }
215 
216  if (phiVal > -LauConstants::pi && phiVal < LauConstants::pi) {
217  phiWithinRange = kTRUE;
218  } else {
219  // Not within the specified range
220  if (phiVal > LauConstants::pi) {
221  phiVal -= LauConstants::twoPi;
222  } else if (phiVal < -LauConstants::pi) {
223  phiVal += LauConstants::twoPi;
224  }
225  }
226  }
227 
228  // A further problem can occur when the generated phase is close to -pi or pi.
229  // The phase can wrap over to the other end of the scale -
230  // this leads to artificially large pulls so we wrap it back.
231  Double_t diff = deltaVal - genDelta;
232  if (diff > LauConstants::pi) {
233  deltaVal -= LauConstants::twoPi;
234  } else if (diff < -LauConstants::pi) {
235  deltaVal += LauConstants::twoPi;
236  }
237 
238  diff = phiVal - genPhi;
239  if (diff > LauConstants::pi) {
240  phiVal -= LauConstants::twoPi;
241  } else if (diff < -LauConstants::pi) {
242  phiVal += LauConstants::twoPi;
243  }
244 
245  // finally store the new values in the parameters
246  // and update the pulls
247  a_->value(aVal); a_->updatePull();
248  b_->value(bVal); b_->updatePull();
249  delta_->value(deltaVal); delta_->updatePull();
250  phi_->value(phiVal); phi_->updatePull();
251 }
252 
254 {
255  LauComplex aTerm(a_->unblindValue()*TMath::Cos(delta_->unblindValue()), a_->unblindValue()*TMath::Sin(delta_->unblindValue()));
256  LauComplex bTerm(b_->unblindValue()*TMath::Cos(phi_->unblindValue()), b_->unblindValue()*TMath::Sin(phi_->unblindValue()));
258  particleCoeff_ += bTerm;
259  particleCoeff_ *= aTerm;
260  return particleCoeff_;
261 }
262 
264 {
265  LauComplex aTerm(a_->unblindValue()*TMath::Cos(delta_->unblindValue()), a_->unblindValue()*TMath::Sin(delta_->unblindValue()));
266  LauComplex bTerm(b_->unblindValue()*TMath::Cos(phi_->unblindValue()), b_->unblindValue()*TMath::Sin(phi_->unblindValue()));
268  antiparticleCoeff_ -= bTerm;
269  antiparticleCoeff_ *= aTerm;
270  return antiparticleCoeff_;
271 }
272 
273 void LauBelleCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
274 {
275  LauComplex sum = coeff + coeffBar;
276  LauComplex diff = coeff - coeffBar;
277  LauComplex ratio = diff / sum;
278 
279  Double_t aVal( 0.5 * sum.abs() );
280  Double_t deltaVal( sum.arg() );
281  Double_t bVal( ratio.abs() );
282  Double_t phiVal( ratio.arg() );
283 
284  a_->value( aVal );
285  delta_->value( deltaVal );
286  b_->value( bVal );
287  phi_->value( phiVal );
288 
289  if ( init ) {
290  a_->genValue( aVal );
291  delta_->genValue( deltaVal );
292  b_->genValue( bVal );
293  phi_->genValue( phiVal );
294 
295  a_->initValue( aVal );
296  delta_->initValue( deltaVal );
297  b_->initValue( bVal );
298  phi_->initValue( phiVal );
299  }
300 }
301 
303 {
304  // set the name
305  TString parName(this->baseName()); parName += "_ACP";
306  acp_.name(parName);
307 
308  // work out the ACP value
309  Double_t value = (-2.0*b_->value()*TMath::Cos(phi_->value()))/(1.0+b_->value()*b_->value());
310 
311  // is it fixed?
312  Bool_t fixed = b_->fixed() && phi_->fixed();
313  acp_.fixed(fixed);
314 
315  // we can't work out the error without the covariance matrix
316  Double_t error(0.0);
317 
318  // set the value and error
319  acp_.valueAndErrors(value,error);
320 
321  return acp_;
322 }
323 
324 LauAbsCoeffSet* LauBelleCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
325 {
326  LauAbsCoeffSet* clone(0);
327  if ( cloneOption == All || cloneOption == TiePhase || cloneOption == TieMagnitude || cloneOption == TieCPPars ) {
328  clone = new LauBelleCPCoeffSet( *this, cloneOption, constFactor );
329  clone->name( newName );
330  } else {
331  std::cerr << "ERROR in LauBelleCPCoeffSet::createClone : Invalid clone option" << std::endl;
332  }
333  return clone;
334 }
335 
virtual void randomiseInitValues()
Randomise the starting values of the parameters for a fit.
static Double_t maxPhase_
Maximum allowed value of phase parameters.
Class for defining a complex coefficient using the Belle CP convention. Holds a set of real values th...
Bool_t fixed() const
Check whether the parameter is fixed or floated.
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.
const Double_t twoPi
Two times Pi.
Definition: LauConstants.hh:93
ClassImp(LauAbsCoeffSet)
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
const TString & name() const
The parameter name.
virtual std::vector< LauParameter * > getParameters()
Retrieve the parameters of the coefficient, e.g. so that they can be loaded into a fit...
virtual LauAbsCoeffSet * createClone(const TString &newName, CloneOption cloneOption=All, Double_t constFactor=1.0)
Create a clone of the coefficient set.
File containing declaration of LauPrint class.
LauBelleCPCoeffSet(const TString &compName, Double_t a, Double_t delta, Double_t b, Double_t phi, Bool_t aFixed, Bool_t deltaFixed, Bool_t bFixed, Bool_t phiFixed, Bool_t bSecondStage=kFALSE, Bool_t phiSecondStage=kFALSE)
Constructor.
Class to define various output print commands.
Definition: LauPrint.hh:29
LauParameter * delta_
The strong phase.
CloneOption
Options for cloning operation.
Bool_t clone() const
Check whether is a clone or not.
Bool_t blind() const
The blinding state.
LauParameter * phi_
The weak phase.
LauParameter * a_
The magnitude a.
static Double_t maxMagnitude_
Maximum allowed value of magnitude parameters.
virtual void printTableHeading(std::ostream &stream) const
Print the column headings for a results table.
File containing declaration of LauParameter class.
LauParameter * b_
The magnitude b.
Bool_t secondStage() const
Check whether the parameter should be floated only in the second stage of a two stage fit...
File containing declaration of LauComplex class.
const TString & blindingString() const
Obtain the blinding string.
Definition: LauBlind.hh:62
virtual LauParameter acp()
Calculate the CP asymmetry.
Double_t error() const
The error on the parameter.
const Double_t pi
Pi.
Definition: LauConstants.hh:89
Class for defining the abstract interface for complex coefficient classes.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:35
void valueAndErrors(Double_t newValue, Double_t newError, Double_t newNegError=0.0, Double_t newPosError=0.0)
Set the value and errors on the parameter.
virtual void printParValues() const
Print the current values of the parameters.
File containing declaration of LauBelleCPCoeffSet class.
const LauBlind * blinder() const
Access the blinder object.
void setRealImagPart(Double_t realpart, Double_t imagpart)
Set both real and imaginary part.
Definition: LauComplex.hh:314
virtual void printTableRow(std::ostream &stream) const
Print the parameters of the complex coefficient as a row in the results table.
virtual void setCoeffValues(const LauComplex &coeff, const LauComplex &coeffBar, Bool_t init)
Set the parameters based on the complex coefficients for particles and antiparticles.
Double_t initValue() const
The initial value of the parameter.
void blindParameter(const TString &blindingString, const Double_t width)
Blind the parameter.
virtual const LauComplex & antiparticleCoeff()
Retrieve the complex coefficient for an antiparticle.
File containing LauConstants namespace.
void printFormat(std::ostream &stream, Double_t value) const
Method to choose the printing format to a specified level of precision.
Definition: LauPrint.cc:32
Double_t unblindValue() const
The unblinded value of the parameter.
Class for defining a complex number.
Definition: LauComplex.hh:47
void updatePull()
Call to update the bias and pull values.
Double_t arg() const
Obtain the phase angle of the complex number.
Definition: LauComplex.hh:241
LauParameter * createClone(Double_t constFactor=1.0)
Method to create a clone from the parent parameter using the copy constructor.
virtual TString name() const
Retrieve the name of the coefficient set.
virtual void finaliseValues()
Make sure values are in &quot;standard&quot; ranges, e.g. phases should be between -pi and pi.
LauComplex particleCoeff_
The particle complex coefficient.
static Double_t minPhase_
Minimum allowed value of phase parameters.
Double_t value() const
The value of the parameter.
static Double_t minMagnitude_
Minimum allowed value of magnitude parameters.
Double_t abs() const
Obtain the absolute value of the complex number.
Definition: LauComplex.hh:223
Double_t blindingWidth() const
Obtain the Gaussian width.
Definition: LauBlind.hh:68
virtual const TString & baseName() const
Retrieve the base name of the coefficient set.
Class for blinding and unblinding a number based on a blinding string.
Definition: LauBlind.hh:28
static TRandom * getRandomiser()
Access the randomiser.
Double_t genValue() const
The value generated for the parameter.
LauParameter acp_
The CP asymmetry.