laura is hosted by Hepforge, IPPP Durham
Laura++  v3r5
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauBelleCPCoeffSet.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2006 University of Warwick
4 
5 Licensed under the Apache License, Version 2.0 (the "License");
6 you may not use this file except in compliance with the License.
7 You may obtain a copy of the License at
8 
9  http://www.apache.org/licenses/LICENSE-2.0
10 
11 Unless required by applicable law or agreed to in writing, software
12 distributed under the License is distributed on an "AS IS" BASIS,
13 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 See the License for the specific language governing permissions and
15 limitations under the License.
16 */
17 
18 /*
19 Laura++ package authors:
20 John Back
21 Paul Harrison
22 Thomas Latham
23 */
24 
29 #include <iostream>
30 #include <fstream>
31 #include <vector>
32 
33 #include "TMath.h"
34 #include "TRandom.h"
35 
36 #include "LauComplex.hh"
37 #include "LauConstants.hh"
38 #include "LauBelleCPCoeffSet.hh"
39 #include "LauParameter.hh"
40 #include "LauPrint.hh"
41 
43 
44 
45 LauBelleCPCoeffSet::LauBelleCPCoeffSet(const TString& compName, Double_t a, Double_t delta, Double_t b, Double_t phi,
46  Bool_t aFixed, Bool_t deltaFixed, Bool_t bFixed, Bool_t phiFixed, Bool_t bSecondStage, Bool_t phiSecondStage) :
47  LauAbsCoeffSet(compName),
48  a_(new LauParameter("A", a, minMagnitude_, maxMagnitude_, aFixed)),
49  b_(new LauParameter("B", b, minMagnitude_, maxMagnitude_, bFixed)),
50  delta_(new LauParameter("Delta", delta, minPhase_, maxPhase_, deltaFixed)),
51  phi_(new LauParameter("Phi", phi, minPhase_, maxPhase_, phiFixed)),
52  particleCoeff_(0.0,0.0),
53  antiparticleCoeff_(0.0,0.0),
54  acp_("ACP", (-2.0*b*TMath::Cos(phi))/(1.0+b*b), -1.0, 1.0, bFixed&&phiFixed)
55 {
56  if (bSecondStage && !bFixed) {
57  b_->secondStage(kTRUE);
58  b_->initValue(0.0);
59  }
60  if (phiSecondStage && !phiFixed) {
61  phi_->secondStage(kTRUE);
62  phi_->initValue(0.0);
63  }
64 }
65 
66 LauBelleCPCoeffSet::LauBelleCPCoeffSet(const LauBelleCPCoeffSet& rhs, CloneOption cloneOption, Double_t constFactor) : LauAbsCoeffSet(rhs.name()),
67  a_(0),
68  b_(0),
69  delta_(0),
70  phi_(0),
71  particleCoeff_( rhs.particleCoeff_ ),
72  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
73  acp_( rhs.acp_ )
74 {
75  if ( cloneOption == All || cloneOption == TieMagnitude ) {
76  a_ = rhs.a_->createClone(constFactor);
77  } else {
78  a_ = new LauParameter("A", rhs.a_->value(), minMagnitude_, maxMagnitude_, rhs.a_->fixed());
79  if ( rhs.a_->blind() ) {
80  const LauBlind* blinder = rhs.a_->blinder();
81  a_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
82  }
83  }
84 
85  if ( cloneOption == All || cloneOption == TieCPPars ) {
86  b_ = rhs.b_->createClone(constFactor);
87  } else {
88  b_ = new LauParameter("B", rhs.b_->value(), minMagnitude_, maxMagnitude_, rhs.b_->fixed());
89  if ( rhs.b_->blind() ) {
90  const LauBlind* blinder = rhs.b_->blinder();
91  b_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
92  }
93  }
94 
95  if ( cloneOption == All || cloneOption == TiePhase ) {
96  delta_ = rhs.delta_->createClone(constFactor);
97  } else {
98  delta_ = new LauParameter("Delta", rhs.delta_->value(), minPhase_, maxPhase_, rhs.delta_->fixed());
99  if ( rhs.delta_->blind() ) {
100  const LauBlind* blinder = rhs.delta_->blinder();
101  delta_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
102  }
103  }
104 
105  if ( cloneOption == All || cloneOption == TieCPPars ) {
106  phi_ = rhs.phi_->createClone(constFactor);
107  } else {
108  phi_ = new LauParameter("Phi", rhs.phi_->value(), minPhase_, maxPhase_, rhs.phi_->fixed());
109  if ( rhs.phi_->blind() ) {
110  const LauBlind* blinder = rhs.phi_->blinder();
111  phi_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
112  }
113  }
114 }
115 
116 std::vector<LauParameter*> LauBelleCPCoeffSet::getParameters()
117 {
118  std::vector<LauParameter*> pars;
119  pars.push_back(a_);
120  pars.push_back(b_);
121  pars.push_back(delta_);
122  pars.push_back(phi_);
123  return pars;
124 }
125 
127 {
128  std::cout<<"INFO in LauBelleCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
129  std::cout<<"a-magnitude = "<<a_->value()<<",\t";
130  std::cout<<"delta = "<<delta_->value()<<",\t";
131  std::cout<<"b-magnitude = "<<b_->value()<<",\t";
132  std::cout<<"phi = "<<phi_->value()<<"."<<std::endl;
133 }
134 
135 void LauBelleCPCoeffSet::printTableHeading(std::ostream& stream) const
136 {
137  stream<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
138  stream<<"\\hline"<<std::endl;
139  stream<<"Component & a-Magnitude & delta & b-Magnitude & phi \\\\"<<std::endl;
140  stream<<"\\hline"<<std::endl;
141 }
142 
143 void LauBelleCPCoeffSet::printTableRow(std::ostream& stream) const
144 {
145  LauPrint print;
146  TString resName = this->name();
147  resName = resName.ReplaceAll("_", "\\_");
148  stream<<resName<<" & $";
149  print.printFormat(stream, a_->value());
150  stream<<" \\pm ";
151  print.printFormat(stream, a_->error());
152  stream<<"$ & $";
153  print.printFormat(stream, delta_->value());
154  stream<<" \\pm ";
155  print.printFormat(stream, delta_->error());
156  stream<<"$ & $";
157  print.printFormat(stream, b_->value());
158  stream<<" \\pm ";
159  print.printFormat(stream, b_->error());
160  stream<<"$ & $";
161  print.printFormat(stream, phi_->value());
162  stream<<" \\pm ";
163  print.printFormat(stream, phi_->error());
164  stream<<"$ \\\\"<<std::endl;
165 }
166 
168 {
169  if (a_->fixed() == kFALSE) {
170  // Choose an a-magnitude between 0.0 and 2.0
171  Double_t mag = LauAbsCoeffSet::getRandomiser()->Rndm()*2.0;
172  a_->initValue(mag); a_->value(mag);
173  }
174  if (b_->fixed() == kFALSE && b_->secondStage() == kFALSE) {
175  // Choose a b-magnitude between 0.0 and 0.1
176  Double_t mag = LauAbsCoeffSet::getRandomiser()->Rndm()*0.1;
177  b_->initValue(mag); b_->value(mag);
178  }
179  if (delta_->fixed() == kFALSE) {
180  // Choose a phase between +- pi
182  delta_->initValue(phase); delta_->value(phase);
183  }
184  if (phi_->fixed() == kFALSE && phi_->secondStage() == kFALSE) {
185  // Choose a phase between +- pi
187  phi_->initValue(phase); phi_->value(phase);
188  }
189 }
190 
192 {
193  // retrieve the current values from the parameters
194  Double_t aVal = a_->value();
195  Double_t bVal = b_->value();
196  Double_t deltaVal = delta_->value();
197  Double_t phiVal = phi_->value();
198  Double_t genDelta = delta_->genValue();
199  Double_t genPhi = phi_->genValue();
200 
201  // Check whether we have a negative "a" magnitude.
202  // If so make it positive and add pi to the "delta" phase.
203  if (aVal < 0.0) {
204  aVal *= -1.0;
205  deltaVal += LauConstants::pi;
206  }
207 
208  // Check whether we have a negative "b" magnitude.
209  // If so make it positive and add pi to the "phi" phase.
210  if (bVal < 0.0) {
211  bVal *= -1.0;
212  phiVal += LauConstants::pi;
213  }
214 
215  // Check now whether the phases lies in the right range (-pi to pi).
216  Bool_t deltaWithinRange(kFALSE);
217  Bool_t phiWithinRange(kFALSE);
218  while (deltaWithinRange == kFALSE && phiWithinRange == kFALSE) {
219  if (deltaVal > -LauConstants::pi && deltaVal < LauConstants::pi) {
220  deltaWithinRange = kTRUE;
221  } else {
222  // Not within the specified range
223  if (deltaVal > LauConstants::pi) {
224  deltaVal -= LauConstants::twoPi;
225  } else if (deltaVal < -LauConstants::pi) {
226  deltaVal += LauConstants::twoPi;
227  }
228  }
229 
230  if (phiVal > -LauConstants::pi && phiVal < LauConstants::pi) {
231  phiWithinRange = kTRUE;
232  } else {
233  // Not within the specified range
234  if (phiVal > LauConstants::pi) {
235  phiVal -= LauConstants::twoPi;
236  } else if (phiVal < -LauConstants::pi) {
237  phiVal += LauConstants::twoPi;
238  }
239  }
240  }
241 
242  // A further problem can occur when the generated phase is close to -pi or pi.
243  // The phase can wrap over to the other end of the scale -
244  // this leads to artificially large pulls so we wrap it back.
245  Double_t diff = deltaVal - genDelta;
246  if (diff > LauConstants::pi) {
247  deltaVal -= LauConstants::twoPi;
248  } else if (diff < -LauConstants::pi) {
249  deltaVal += LauConstants::twoPi;
250  }
251 
252  diff = phiVal - genPhi;
253  if (diff > LauConstants::pi) {
254  phiVal -= LauConstants::twoPi;
255  } else if (diff < -LauConstants::pi) {
256  phiVal += LauConstants::twoPi;
257  }
258 
259  // finally store the new values in the parameters
260  // and update the pulls
261  a_->value(aVal); a_->updatePull();
262  b_->value(bVal); b_->updatePull();
263  delta_->value(deltaVal); delta_->updatePull();
264  phi_->value(phiVal); phi_->updatePull();
265 }
266 
268 {
269  LauComplex aTerm(a_->unblindValue()*TMath::Cos(delta_->unblindValue()), a_->unblindValue()*TMath::Sin(delta_->unblindValue()));
270  LauComplex bTerm(b_->unblindValue()*TMath::Cos(phi_->unblindValue()), b_->unblindValue()*TMath::Sin(phi_->unblindValue()));
272  particleCoeff_ += bTerm;
273  particleCoeff_ *= aTerm;
274  return particleCoeff_;
275 }
276 
278 {
279  LauComplex aTerm(a_->unblindValue()*TMath::Cos(delta_->unblindValue()), a_->unblindValue()*TMath::Sin(delta_->unblindValue()));
280  LauComplex bTerm(b_->unblindValue()*TMath::Cos(phi_->unblindValue()), b_->unblindValue()*TMath::Sin(phi_->unblindValue()));
282  antiparticleCoeff_ -= bTerm;
283  antiparticleCoeff_ *= aTerm;
284  return antiparticleCoeff_;
285 }
286 
287 void LauBelleCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
288 {
289  LauComplex sum = coeff + coeffBar;
290  LauComplex diff = coeff - coeffBar;
291  LauComplex ratio = diff / sum;
292 
293  Double_t aVal( 0.5 * sum.abs() );
294  Double_t deltaVal( sum.arg() );
295  Double_t bVal( ratio.abs() );
296  Double_t phiVal( ratio.arg() );
297 
298  a_->value( aVal );
299  delta_->value( deltaVal );
300  b_->value( bVal );
301  phi_->value( phiVal );
302 
303  if ( init ) {
304  a_->genValue( aVal );
305  delta_->genValue( deltaVal );
306  b_->genValue( bVal );
307  phi_->genValue( phiVal );
308 
309  a_->initValue( aVal );
310  delta_->initValue( deltaVal );
311  b_->initValue( bVal );
312  phi_->initValue( phiVal );
313  }
314 }
315 
317 {
318  // set the name
319  TString parName(this->baseName()); parName += "_ACP";
320  acp_.name(parName);
321 
322  // work out the ACP value
323  Double_t value = (-2.0*b_->value()*TMath::Cos(phi_->value()))/(1.0+b_->value()*b_->value());
324 
325  // is it fixed?
326  Bool_t fixed = b_->fixed() && phi_->fixed();
327  acp_.fixed(fixed);
328 
329  // we can't work out the error without the covariance matrix
330  Double_t error(0.0);
331 
332  // set the value and error
333  acp_.valueAndErrors(value,error);
334 
335  return acp_;
336 }
337 
338 LauAbsCoeffSet* LauBelleCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
339 {
340  LauAbsCoeffSet* clone(0);
341  if ( cloneOption == All || cloneOption == TiePhase || cloneOption == TieMagnitude || cloneOption == TieCPPars ) {
342  clone = new LauBelleCPCoeffSet( *this, cloneOption, constFactor );
343  clone->name( newName );
344  } else {
345  std::cerr << "ERROR in LauBelleCPCoeffSet::createClone : Invalid clone option" << std::endl;
346  }
347  return clone;
348 }
349 
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.
ClassImp(LauAbsCoeffSet)
LauParameter()
Default constructor.
Definition: LauParameter.cc:44
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:43
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:76
virtual LauParameter acp()
Calculate the CP asymmetry.
Double_t error() const
The error on the parameter.
const Double_t pi
Pi.
Class for defining the abstract interface for complex coefficient classes.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:49
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:328
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:46
Double_t unblindValue() const
The unblinded value of the parameter.
Class for defining a complex number.
Definition: LauComplex.hh:61
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:255
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:237
Double_t blindingWidth() const
Obtain the Gaussian width.
Definition: LauBlind.hh:82
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:42
static TRandom * getRandomiser()
Access the randomiser.
Double_t genValue() const
The value generated for the parameter.
LauParameter acp_
The CP asymmetry.