laura is hosted by Hepforge, IPPP Durham
Laura++  v3r4
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauCleoCPCoeffSet.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 "LauCleoCPCoeffSet.hh"
37 #include "LauComplex.hh"
38 #include "LauConstants.hh"
39 #include "LauParameter.hh"
40 #include "LauPrint.hh"
41 
43 
44 
45 LauCleoCPCoeffSet::LauCleoCPCoeffSet(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_( (a+b)*TMath::Cos(delta+phi), (a+b)*TMath::Sin(delta+phi) ),
53  antiparticleCoeff_( (a-b)*TMath::Cos(delta-phi), (a-b)*TMath::Sin(delta-phi) ),
54  acp_("ACP", (-2.0*a*b)/(a*a+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 LauCleoCPCoeffSet::LauCleoCPCoeffSet(const LauCleoCPCoeffSet& 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*> LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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  bVal *= -1.0;
206  deltaVal += LauConstants::pi;
207  }
208 
209  // Check now whether the phases lies in the right range (-pi to pi).
210  Bool_t deltaWithinRange(kFALSE);
211  Bool_t phiWithinRange(kFALSE);
212  while (deltaWithinRange == kFALSE && phiWithinRange == kFALSE) {
213  if (deltaVal > -LauConstants::pi && deltaVal < LauConstants::pi) {
214  deltaWithinRange = kTRUE;
215  } else {
216  // Not within the specified range
217  if (deltaVal > LauConstants::pi) {
218  deltaVal -= LauConstants::twoPi;
219  } else if (deltaVal < -LauConstants::pi) {
220  deltaVal += LauConstants::twoPi;
221  }
222  }
223 
224  if (phiVal > -LauConstants::pi && phiVal < LauConstants::pi) {
225  phiWithinRange = kTRUE;
226  } else {
227  // Not within the specified range
228  if (phiVal > LauConstants::pi) {
229  phiVal -= LauConstants::twoPi;
230  } else if (phiVal < -LauConstants::pi) {
231  phiVal += LauConstants::twoPi;
232  }
233  }
234  }
235 
236  // A further problem can occur when the generated phase is close to -pi or pi.
237  // The phase can wrap over to the other end of the scale -
238  // this leads to artificially large pulls so we wrap it back.
239  Double_t diff = deltaVal - genDelta;
240  if (diff > LauConstants::pi) {
241  deltaVal -= LauConstants::twoPi;
242  } else if (diff < -LauConstants::pi) {
243  deltaVal += LauConstants::twoPi;
244  }
245 
246  diff = phiVal - genPhi;
247  if (diff > LauConstants::pi) {
248  phiVal -= LauConstants::twoPi;
249  } else if (diff < -LauConstants::pi) {
250  phiVal += LauConstants::twoPi;
251  }
252 
253  // finally store the new values in the parameters
254  // and update the pulls
255  a_->value(aVal); a_->updatePull();
256  b_->value(bVal); b_->updatePull();
257  delta_->value(deltaVal); delta_->updatePull();
258  phi_->value(phiVal); phi_->updatePull();
259 }
260 
262 {
263  Double_t magnitude = a_->unblindValue() + b_->unblindValue();
264  Double_t phase = delta_->unblindValue() + phi_->unblindValue();
265  particleCoeff_.setRealImagPart(magnitude*TMath::Cos(phase), magnitude*TMath::Sin(phase));
266  return particleCoeff_;
267 }
268 
270 {
271  Double_t magnitude = a_->unblindValue() - b_->unblindValue();
272  Double_t phase = delta_->unblindValue() - phi_->unblindValue();
273  antiparticleCoeff_.setRealImagPart(magnitude*TMath::Cos(phase), magnitude*TMath::Sin(phase));
274  return antiparticleCoeff_;
275 }
276 
277 void LauCleoCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
278 {
279  Double_t mag = coeff.abs();
280  Double_t magBar = coeffBar.abs();
281  Double_t phase = coeff.arg();
282  Double_t phaseBar = coeffBar.arg();
283 
284  Double_t aVal( 0.5 * ( mag + magBar ) );
285  Double_t deltaVal( 0.5 * ( phase + phaseBar ) );
286  Double_t bVal( 0.5 * ( mag - magBar ) );
287  Double_t phiVal( 0.5 * ( phase - phaseBar ) );
288 
289  a_->value( aVal );
290  delta_->value( deltaVal );
291  b_->value( bVal );
292  phi_->value( phiVal );
293 
294  if ( init ) {
295  a_->genValue( aVal );
296  delta_->genValue( deltaVal );
297  b_->genValue( bVal );
298  phi_->genValue( phiVal );
299 
300  a_->initValue( aVal );
301  delta_->initValue( deltaVal );
302  b_->initValue( bVal );
303  phi_->initValue( phiVal );
304  }
305 }
306 
308 {
309  // set the name
310  TString parName(this->baseName()); parName += "_ACP";
311  acp_.name(parName);
312 
313  // work out the ACP value
314  Double_t numer = -2.0*a_->value()*b_->value();
315  Double_t denom = a_->value()*a_->value()+b_->value()*b_->value();
316  Double_t value = numer/denom;
317 
318  // is it fixed?
319  Bool_t fixed = a_->fixed() && b_->fixed();
320  acp_.fixed(fixed);
321 
322  // we can't work out the error without the covariance matrix
323  Double_t error(0.0);
324 
325  // set the value and error
326  acp_.valueAndErrors(value,error);
327 
328  return acp_;
329 }
330 
331 LauAbsCoeffSet* LauCleoCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
332 {
333  LauAbsCoeffSet* clone(0);
334  if ( cloneOption == All || cloneOption == TiePhase || cloneOption == TieMagnitude || cloneOption == TieCPPars ) {
335  clone = new LauCleoCPCoeffSet( *this, cloneOption, constFactor );
336  clone->name( newName );
337  } else {
338  std::cerr << "ERROR in LauCleoCPCoeffSet::createClone : Invalid clone option" << std::endl;
339  }
340  return clone;
341 }
342 
virtual void printParValues() const
Print the current values of the parameters.
LauParameter acp_
The CP asymmetry.
static Double_t maxPhase_
Maximum allowed value of phase parameters.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
File containing declaration of LauCleoCPCoeffSet class.
const Double_t twoPi
Two times Pi.
ClassImp(LauAbsCoeffSet)
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.
LauParameter()
Default constructor.
Definition: LauParameter.cc:44
const TString & name() const
The parameter name.
LauParameter * delta_
The strong phase.
File containing declaration of LauPrint class.
Class to define various output print commands.
Definition: LauPrint.hh:43
LauCleoCPCoeffSet(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 for defining a complex coefficient using the Cleo CP convention.
LauComplex particleCoeff_
The particle complex coefficient.
LauParameter * phi_
The weak phase.
CloneOption
Options for cloning operation.
virtual const LauComplex & antiparticleCoeff()
Retrieve the complex coefficient for an antiparticle.
Bool_t clone() const
Check whether is a clone or not.
virtual LauParameter acp()
Calculate the CP asymmetry.
Bool_t blind() const
The blinding state.
virtual void printTableHeading(std::ostream &stream) const
Print the column headings for a results table.
static Double_t maxMagnitude_
Maximum allowed value of magnitude parameters.
File containing declaration of LauParameter class.
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
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.
LauParameter * b_
The magnitude b.
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 LauAbsCoeffSet * createClone(const TString &newName, CloneOption cloneOption=All, Double_t constFactor=1.0)
Create a clone of the coefficient set.
Double_t initValue() const
The initial value of the parameter.
void blindParameter(const TString &blindingString, const Double_t width)
Blind the parameter.
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.
virtual void randomiseInitValues()
Randomise the starting values of the parameters for a fit.
virtual void setCoeffValues(const LauComplex &coeff, const LauComplex &coeffBar, Bool_t init)
Set the parameters based on the complex coefficients for particles and antiparticles.
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.
static Double_t minPhase_
Minimum allowed value of phase parameters.
Double_t value() const
The value of the parameter.
virtual void printTableRow(std::ostream &stream) const
Print the parameters of the complex coefficient as a row in the results table.
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
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
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.
virtual void finaliseValues()
Make sure values are in &quot;standard&quot; ranges, e.g. phases should be between -pi and pi.
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.
virtual std::vector< LauParameter * > getParameters()
Retrieve the parameters of the coefficient, e.g. so that they can be loaded into a fit...
LauParameter * a_
The magnitude a.