laura is hosted by Hepforge, IPPP Durham
Laura++  v3r2
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauCleoCPCoeffSet.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 "LauCleoCPCoeffSet.hh"
23 #include "LauComplex.hh"
24 #include "LauConstants.hh"
25 #include "LauParameter.hh"
26 #include "LauPrint.hh"
27 
29 
30 
31 LauCleoCPCoeffSet::LauCleoCPCoeffSet(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_( (a+b)*TMath::Cos(delta+phi), (a+b)*TMath::Sin(delta+phi) ),
39  antiparticleCoeff_( (a-b)*TMath::Cos(delta-phi), (a-b)*TMath::Sin(delta-phi) ),
40  acp_("ACP", (-2.0*a*b)/(a*a+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 LauCleoCPCoeffSet::LauCleoCPCoeffSet(const LauCleoCPCoeffSet& 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*> LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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  bVal *= -1.0;
192  deltaVal += LauConstants::pi;
193  }
194 
195  // Check now whether the phases lies in the right range (-pi to pi).
196  Bool_t deltaWithinRange(kFALSE);
197  Bool_t phiWithinRange(kFALSE);
198  while (deltaWithinRange == kFALSE && phiWithinRange == kFALSE) {
199  if (deltaVal > -LauConstants::pi && deltaVal < LauConstants::pi) {
200  deltaWithinRange = kTRUE;
201  } else {
202  // Not within the specified range
203  if (deltaVal > LauConstants::pi) {
204  deltaVal -= LauConstants::twoPi;
205  } else if (deltaVal < -LauConstants::pi) {
206  deltaVal += LauConstants::twoPi;
207  }
208  }
209 
210  if (phiVal > -LauConstants::pi && phiVal < LauConstants::pi) {
211  phiWithinRange = kTRUE;
212  } else {
213  // Not within the specified range
214  if (phiVal > LauConstants::pi) {
215  phiVal -= LauConstants::twoPi;
216  } else if (phiVal < -LauConstants::pi) {
217  phiVal += LauConstants::twoPi;
218  }
219  }
220  }
221 
222  // A further problem can occur when the generated phase is close to -pi or pi.
223  // The phase can wrap over to the other end of the scale -
224  // this leads to artificially large pulls so we wrap it back.
225  Double_t diff = deltaVal - genDelta;
226  if (diff > LauConstants::pi) {
227  deltaVal -= LauConstants::twoPi;
228  } else if (diff < -LauConstants::pi) {
229  deltaVal += LauConstants::twoPi;
230  }
231 
232  diff = phiVal - genPhi;
233  if (diff > LauConstants::pi) {
234  phiVal -= LauConstants::twoPi;
235  } else if (diff < -LauConstants::pi) {
236  phiVal += LauConstants::twoPi;
237  }
238 
239  // finally store the new values in the parameters
240  // and update the pulls
241  a_->value(aVal); a_->updatePull();
242  b_->value(bVal); b_->updatePull();
243  delta_->value(deltaVal); delta_->updatePull();
244  phi_->value(phiVal); phi_->updatePull();
245 }
246 
248 {
249  Double_t magnitude = a_->unblindValue() + b_->unblindValue();
250  Double_t phase = delta_->unblindValue() + phi_->unblindValue();
251  particleCoeff_.setRealImagPart(magnitude*TMath::Cos(phase), magnitude*TMath::Sin(phase));
252  return particleCoeff_;
253 }
254 
256 {
257  Double_t magnitude = a_->unblindValue() - b_->unblindValue();
258  Double_t phase = delta_->unblindValue() - phi_->unblindValue();
259  antiparticleCoeff_.setRealImagPart(magnitude*TMath::Cos(phase), magnitude*TMath::Sin(phase));
260  return antiparticleCoeff_;
261 }
262 
263 void LauCleoCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
264 {
265  Double_t mag = coeff.abs();
266  Double_t magBar = coeffBar.abs();
267  Double_t phase = coeff.arg();
268  Double_t phaseBar = coeffBar.arg();
269 
270  Double_t aVal( 0.5 * ( mag + magBar ) );
271  Double_t deltaVal( 0.5 * ( phase + phaseBar ) );
272  Double_t bVal( 0.5 * ( mag - magBar ) );
273  Double_t phiVal( 0.5 * ( phase - phaseBar ) );
274 
275  a_->value( aVal );
276  delta_->value( deltaVal );
277  b_->value( bVal );
278  phi_->value( phiVal );
279 
280  if ( init ) {
281  a_->genValue( aVal );
282  delta_->genValue( deltaVal );
283  b_->genValue( bVal );
284  phi_->genValue( phiVal );
285 
286  a_->initValue( aVal );
287  delta_->initValue( deltaVal );
288  b_->initValue( bVal );
289  phi_->initValue( phiVal );
290  }
291 }
292 
294 {
295  // set the name
296  TString parName(this->baseName()); parName += "_ACP";
297  acp_.name(parName);
298 
299  // work out the ACP value
300  Double_t numer = -2.0*a_->value()*b_->value();
301  Double_t denom = a_->value()*a_->value()+b_->value()*b_->value();
302  Double_t value = numer/denom;
303 
304  // is it fixed?
305  Bool_t fixed = a_->fixed() && b_->fixed();
306  acp_.fixed(fixed);
307 
308  // we can't work out the error without the covariance matrix
309  Double_t error(0.0);
310 
311  // set the value and error
312  acp_.valueAndErrors(value,error);
313 
314  return acp_;
315 }
316 
317 LauAbsCoeffSet* LauCleoCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
318 {
319  LauAbsCoeffSet* clone(0);
320  if ( cloneOption == All || cloneOption == TiePhase || cloneOption == TieMagnitude || cloneOption == TieCPPars ) {
321  clone = new LauCleoCPCoeffSet( *this, cloneOption, constFactor );
322  clone->name( newName );
323  } else {
324  std::cerr << "ERROR in LauCleoCPCoeffSet::createClone : Invalid clone option" << std::endl;
325  }
326  return clone;
327 }
328 
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.
Definition: LauConstants.hh:93
ClassImp(LauAbsCoeffSet)
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
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:29
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:62
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.
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:314
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:32
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: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.
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:223
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
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.
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:28
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.