laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0p1
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 #include "LauRandom.hh"
28 
30 
31 
32 LauCleoCPCoeffSet::LauCleoCPCoeffSet(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_( (a+b)*TMath::Cos(delta+phi), (a+b)*TMath::Sin(delta+phi) ),
40  antiparticleCoeff_( (a-b)*TMath::Cos(delta-phi), (a-b)*TMath::Sin(delta-phi) ),
41  acp_("ACP", (-2.0*a*b)/(a*a+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 LauCleoCPCoeffSet::LauCleoCPCoeffSet(const LauCleoCPCoeffSet& 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*> LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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 LauCleoCPCoeffSet::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  bVal *= -1.0;
177  deltaVal += LauConstants::pi;
178  }
179 
180  // Check now whether the phases lies in the right range (-pi to pi).
181  Bool_t deltaWithinRange(kFALSE);
182  Bool_t phiWithinRange(kFALSE);
183  while (deltaWithinRange == kFALSE && phiWithinRange == kFALSE) {
184  if (deltaVal > -LauConstants::pi && deltaVal < LauConstants::pi) {
185  deltaWithinRange = kTRUE;
186  } else {
187  // Not within the specified range
188  if (deltaVal > LauConstants::pi) {
189  deltaVal -= LauConstants::twoPi;
190  } else if (deltaVal < -LauConstants::pi) {
191  deltaVal += LauConstants::twoPi;
192  }
193  }
194 
195  if (phiVal > -LauConstants::pi && phiVal < LauConstants::pi) {
196  phiWithinRange = kTRUE;
197  } else {
198  // Not within the specified range
199  if (phiVal > LauConstants::pi) {
200  phiVal -= LauConstants::twoPi;
201  } else if (phiVal < -LauConstants::pi) {
202  phiVal += LauConstants::twoPi;
203  }
204  }
205  }
206 
207  // A further problem can occur when the generated phase is close to -pi or pi.
208  // The phase can wrap over to the other end of the scale -
209  // this leads to artificially large pulls so we wrap it back.
210  Double_t diff = deltaVal - genDelta;
211  if (diff > LauConstants::pi) {
212  deltaVal -= LauConstants::twoPi;
213  } else if (diff < -LauConstants::pi) {
214  deltaVal += LauConstants::twoPi;
215  }
216 
217  diff = phiVal - genPhi;
218  if (diff > LauConstants::pi) {
219  phiVal -= LauConstants::twoPi;
220  } else if (diff < -LauConstants::pi) {
221  phiVal += LauConstants::twoPi;
222  }
223 
224  // finally store the new values in the parameters
225  // and update the pulls
226  a_->value(aVal); a_->updatePull();
227  b_->value(bVal); b_->updatePull();
228  delta_->value(deltaVal); delta_->updatePull();
229  phi_->value(phiVal); phi_->updatePull();
230 }
231 
233 {
234  Double_t magnitude = a_->value() + b_->value();
235  Double_t phase = delta_->value() + phi_->value();
236  particleCoeff_.setRealImagPart(magnitude*TMath::Cos(phase), magnitude*TMath::Sin(phase));
237  return particleCoeff_;
238 }
239 
241 {
242  Double_t magnitude = a_->value() - b_->value();
243  Double_t phase = delta_->value() - phi_->value();
244  antiparticleCoeff_.setRealImagPart(magnitude*TMath::Cos(phase), magnitude*TMath::Sin(phase));
245  return antiparticleCoeff_;
246 }
247 
248 void LauCleoCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
249 {
250  Double_t mag = coeff.abs();
251  Double_t magBar = coeffBar.abs();
252  Double_t phase = coeff.arg();
253  Double_t phaseBar = coeffBar.arg();
254 
255  Double_t aVal( 0.5 * ( mag + magBar ) );
256  Double_t deltaVal( 0.5 * ( phase + phaseBar ) );
257  Double_t bVal( 0.5 * ( mag - magBar ) );
258  Double_t phiVal( 0.5 * ( phase - phaseBar ) );
259 
260  a_->value( aVal );
261  delta_->value( deltaVal );
262  b_->value( bVal );
263  phi_->value( phiVal );
264 
265  if ( init ) {
266  a_->genValue( aVal );
267  delta_->genValue( deltaVal );
268  b_->genValue( bVal );
269  phi_->genValue( phiVal );
270 
271  a_->initValue( aVal );
272  delta_->initValue( deltaVal );
273  b_->initValue( bVal );
274  phi_->initValue( phiVal );
275  }
276 }
277 
279 {
280  // set the name
281  TString parName(this->baseName()); parName += "_ACP";
282  acp_.name(parName);
283 
284  // work out the ACP value
285  Double_t numer = -2.0*a_->value()*b_->value();
286  Double_t denom = a_->value()*a_->value()+b_->value()*b_->value();
287  Double_t value = numer/denom;
288 
289  // is it fixed?
290  Bool_t fixed = a_->fixed() && b_->fixed();
291  acp_.fixed(fixed);
292 
293  // we can't work out the error without the covariance matrix
294  Double_t error(0.0);
295 
296  // set the value and error
297  acp_.valueAndErrors(value,error);
298 
299  return acp_;
300 }
301 
302 LauAbsCoeffSet* LauCleoCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
303 {
304  LauAbsCoeffSet* clone(0);
305  if ( cloneOption == All || cloneOption == TiePhase || cloneOption == TieMagnitude || cloneOption == TieCPPars ) {
306  clone = new LauCleoCPCoeffSet( *this, cloneOption, constFactor );
307  clone->name( newName );
308  } else {
309  std::cerr << "ERROR in LauCleoCPCoeffSet::createClone : Invalid clone option" << std::endl;
310  }
311  return clone;
312 }
313 
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.
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)
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.
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.
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.
LauParameter * b_
The magnitude b.
File containing LauRandom namespace.
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.
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
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.
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.
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.