laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0
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 #include "LauRandom.hh"
28 
30 
31 
32 LauBelleCPCoeffSet::LauBelleCPCoeffSet(const TString& compName, Double_t a, Double_t delta, Double_t b, Double_t phi,
33  Bool_t aFixed, Bool_t deltaFixed, Bool_t bFixed, Bool_t phiFixed, Bool_t bSecondStage, Bool_t phiSecondStage) :
34  LauAbsCoeffSet(compName),
35  a_(new LauParameter("A", a, minMagnitude_, maxMagnitude_, aFixed)),
36  b_(new LauParameter("B", b, minMagnitude_, maxMagnitude_, bFixed)),
37  delta_(new LauParameter("Delta", delta, minPhase_, maxPhase_, deltaFixed)),
38  phi_(new LauParameter("Phi", phi, minPhase_, maxPhase_, phiFixed)),
39  particleCoeff_(0.0,0.0),
40  antiparticleCoeff_(0.0,0.0),
41  acp_("ACP", (-2.0*b*TMath::Cos(phi))/(1.0+b*b), -1.0, 1.0, bFixed&&phiFixed)
42 {
43  if (bSecondStage && !bFixed) {
44  b_->secondStage(kTRUE);
45  b_->initValue(0.0);
46  }
47  if (phiSecondStage && !phiFixed) {
48  phi_->secondStage(kTRUE);
49  phi_->initValue(0.0);
50  }
51 }
52 
53 LauBelleCPCoeffSet::LauBelleCPCoeffSet(const LauBelleCPCoeffSet& rhs, CloneOption cloneOption, Double_t constFactor) : LauAbsCoeffSet(rhs.name()),
54  a_(0),
55  b_(0),
56  delta_(0),
57  phi_(0),
58  particleCoeff_( rhs.particleCoeff_ ),
59  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
60  acp_( rhs.acp_ )
61 {
62  if ( cloneOption == All || cloneOption == TieMagnitude ) {
63  a_ = rhs.a_->createClone(constFactor);
64  } else {
65  a_ = new LauParameter("A", rhs.a_->value(), minMagnitude_, maxMagnitude_, rhs.a_->fixed());
66  }
67 
68  if ( cloneOption == All || cloneOption == TieCPPars ) {
69  b_ = rhs.b_->createClone(constFactor);
70  } else {
71  b_ = new LauParameter("B", rhs.b_->value(), minMagnitude_, maxMagnitude_, rhs.b_->fixed());
72  }
73 
74  if ( cloneOption == All || cloneOption == TiePhase ) {
75  delta_ = rhs.delta_->createClone(constFactor);
76  } else {
77  delta_ = new LauParameter("Delta", rhs.delta_->value(), minPhase_, maxPhase_, rhs.delta_->fixed());
78  }
79 
80  if ( cloneOption == All || cloneOption == TieCPPars ) {
81  phi_ = rhs.phi_->createClone(constFactor);
82  } else {
83  phi_ = new LauParameter("Phi", rhs.phi_->value(), minPhase_, maxPhase_, rhs.phi_->fixed());
84  }
85 }
86 
87 std::vector<LauParameter*> LauBelleCPCoeffSet::getParameters()
88 {
89  std::vector<LauParameter*> pars;
90  pars.push_back(a_);
91  pars.push_back(b_);
92  pars.push_back(delta_);
93  pars.push_back(phi_);
94  return pars;
95 }
96 
98 {
99  std::cout<<"INFO in LauBelleCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
100  std::cout<<"a-magnitude = "<<a_->value()<<",\t";
101  std::cout<<"delta = "<<delta_->value()<<",\t";
102  std::cout<<"b-magnitude = "<<b_->value()<<",\t";
103  std::cout<<"phi = "<<phi_->value()<<"."<<std::endl;
104 }
105 
106 void LauBelleCPCoeffSet::printTableHeading(std::ostream& stream) const
107 {
108  stream<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
109  stream<<"\\hline"<<std::endl;
110  stream<<"Component & a-Magnitude & delta & b-Magnitude & phi \\\\"<<std::endl;
111  stream<<"\\hline"<<std::endl;
112 }
113 
114 void LauBelleCPCoeffSet::printTableRow(std::ostream& stream) const
115 {
116  LauPrint print;
117  TString resName = this->name();
118  resName = resName.ReplaceAll("_", "\\_");
119  stream<<resName<<" & $";
120  print.printFormat(stream, a_->value());
121  stream<<" \\pm ";
122  print.printFormat(stream, a_->error());
123  stream<<"$ & $";
124  print.printFormat(stream, delta_->value());
125  stream<<" \\pm ";
126  print.printFormat(stream, delta_->error());
127  stream<<"$ & $";
128  print.printFormat(stream, b_->value());
129  stream<<" \\pm ";
130  print.printFormat(stream, b_->error());
131  stream<<"$ & $";
132  print.printFormat(stream, phi_->value());
133  stream<<" \\pm ";
134  print.printFormat(stream, phi_->error());
135  stream<<"$ \\\\"<<std::endl;
136 }
137 
139 {
140  if (a_->fixed() == kFALSE) {
141  // Choose an a-magnitude between 0.0 and 2.0
142  Double_t mag = LauRandom::zeroSeedRandom()->Rndm()*2.0;
143  a_->initValue(mag); a_->value(mag);
144  }
145  if (b_->fixed() == kFALSE && b_->secondStage() == kFALSE) {
146  // Choose a b-magnitude between 0.0 and 0.1
147  Double_t mag = LauRandom::zeroSeedRandom()->Rndm()*0.1;
148  b_->initValue(mag); b_->value(mag);
149  }
150  if (delta_->fixed() == kFALSE) {
151  // Choose a phase between +- pi
153  delta_->initValue(phase); delta_->value(phase);
154  }
155  if (phi_->fixed() == kFALSE && phi_->secondStage() == kFALSE) {
156  // Choose a phase between +- pi
158  phi_->initValue(phase); phi_->value(phase);
159  }
160 }
161 
163 {
164  // retrieve the current values from the parameters
165  Double_t aVal = a_->value();
166  Double_t bVal = b_->value();
167  Double_t deltaVal = delta_->value();
168  Double_t phiVal = phi_->value();
169  Double_t genDelta = delta_->genValue();
170  Double_t genPhi = phi_->genValue();
171 
172  // Check whether we have a negative "a" magnitude.
173  // If so make it positive and add pi to the "delta" phase.
174  if (aVal < 0.0) {
175  aVal *= -1.0;
176  deltaVal += LauConstants::pi;
177  }
178 
179  // Check whether we have a negative "b" magnitude.
180  // If so make it positive and add pi to the "phi" phase.
181  if (bVal < 0.0) {
182  bVal *= -1.0;
183  phiVal += LauConstants::pi;
184  }
185 
186  // Check now whether the phases lies in the right range (-pi to pi).
187  Bool_t deltaWithinRange(kFALSE);
188  Bool_t phiWithinRange(kFALSE);
189  while (deltaWithinRange == kFALSE && phiWithinRange == kFALSE) {
190  if (deltaVal > -LauConstants::pi && deltaVal < LauConstants::pi) {
191  deltaWithinRange = kTRUE;
192  } else {
193  // Not within the specified range
194  if (deltaVal > LauConstants::pi) {
195  deltaVal -= LauConstants::twoPi;
196  } else if (deltaVal < -LauConstants::pi) {
197  deltaVal += LauConstants::twoPi;
198  }
199  }
200 
201  if (phiVal > -LauConstants::pi && phiVal < LauConstants::pi) {
202  phiWithinRange = kTRUE;
203  } else {
204  // Not within the specified range
205  if (phiVal > LauConstants::pi) {
206  phiVal -= LauConstants::twoPi;
207  } else if (phiVal < -LauConstants::pi) {
208  phiVal += LauConstants::twoPi;
209  }
210  }
211  }
212 
213  // A further problem can occur when the generated phase is close to -pi or pi.
214  // The phase can wrap over to the other end of the scale -
215  // this leads to artificially large pulls so we wrap it back.
216  Double_t diff = deltaVal - genDelta;
217  if (diff > LauConstants::pi) {
218  deltaVal -= LauConstants::twoPi;
219  } else if (diff < -LauConstants::pi) {
220  deltaVal += LauConstants::twoPi;
221  }
222 
223  diff = phiVal - genPhi;
224  if (diff > LauConstants::pi) {
225  phiVal -= LauConstants::twoPi;
226  } else if (diff < -LauConstants::pi) {
227  phiVal += LauConstants::twoPi;
228  }
229 
230  // finally store the new values in the parameters
231  // and update the pulls
232  a_->value(aVal); a_->updatePull();
233  b_->value(bVal); b_->updatePull();
234  delta_->value(deltaVal); delta_->updatePull();
235  phi_->value(phiVal); phi_->updatePull();
236 }
237 
239 {
240  LauComplex aTerm(a_->value()*TMath::Cos(delta_->value()), a_->value()*TMath::Sin(delta_->value()));
241  LauComplex bTerm(b_->value()*TMath::Cos(phi_->value()), b_->value()*TMath::Sin(phi_->value()));
243  particleCoeff_ += bTerm;
244  particleCoeff_ *= aTerm;
245  return particleCoeff_;
246 }
247 
249 {
250  LauComplex aTerm(a_->value()*TMath::Cos(delta_->value()), a_->value()*TMath::Sin(delta_->value()));
251  LauComplex bTerm(b_->value()*TMath::Cos(phi_->value()), b_->value()*TMath::Sin(phi_->value()));
253  antiparticleCoeff_ -= bTerm;
254  antiparticleCoeff_ *= aTerm;
255  return antiparticleCoeff_;
256 }
257 
258 void LauBelleCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
259 {
260  LauComplex sum = coeff + coeffBar;
261  LauComplex diff = coeff - coeffBar;
262  LauComplex ratio = diff / sum;
263 
264  Double_t aVal( 0.5 * sum.abs() );
265  Double_t deltaVal( sum.arg() );
266  Double_t bVal( ratio.abs() );
267  Double_t phiVal( ratio.arg() );
268 
269  a_->value( aVal );
270  delta_->value( deltaVal );
271  b_->value( bVal );
272  phi_->value( phiVal );
273 
274  if ( init ) {
275  a_->genValue( aVal );
276  delta_->genValue( deltaVal );
277  b_->genValue( bVal );
278  phi_->genValue( phiVal );
279 
280  a_->initValue( aVal );
281  delta_->initValue( deltaVal );
282  b_->initValue( bVal );
283  phi_->initValue( phiVal );
284  }
285 }
286 
288 {
289  // set the name
290  TString parName(this->baseName()); parName += "_ACP";
291  acp_.name(parName);
292 
293  // work out the ACP value
294  Double_t value = (-2.0*b_->value()*TMath::Cos(phi_->value()))/(1.0+b_->value()*b_->value());
295 
296  // is it fixed?
297  Bool_t fixed = b_->fixed() && phi_->fixed();
298  acp_.fixed(fixed);
299 
300  // we can't work out the error without the covariance matrix
301  Double_t error(0.0);
302 
303  // set the value and error
304  acp_.valueAndErrors(value,error);
305 
306  return acp_;
307 }
308 
309 LauAbsCoeffSet* LauBelleCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
310 {
311  LauAbsCoeffSet* clone(0);
312  if ( cloneOption == All || cloneOption == TiePhase || cloneOption == TieMagnitude || cloneOption == TieCPPars ) {
313  clone = new LauBelleCPCoeffSet( *this, cloneOption, constFactor );
314  clone->name( newName );
315  } else {
316  std::cerr << "ERROR in LauBelleCPCoeffSet::createClone : Invalid clone option" << std::endl;
317  }
318  return clone;
319 }
320 
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.
TRandom * zeroSeedRandom()
Access the singleton random number generator with seed set from machine clock time (within +-1 sec)...
Definition: LauRandom.cc:30
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.
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.
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:34
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 LauRandom namespace.
File containing declaration of LauBelleCPCoeffSet class.
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.
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
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
virtual const TString & baseName() const
Retrieve the base name of the coefficient set.
Double_t genValue() const
The value generated for the parameter.
LauParameter acp_
The CP asymmetry.