laura is hosted by Hepforge, IPPP Durham
Laura++  v3r1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauCartesianCPCoeffSet.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 
23 #include "LauComplex.hh"
24 #include "LauConstants.hh"
25 #include "LauParameter.hh"
26 #include "LauPrint.hh"
27 #include "LauRandom.hh"
28 
30 
31 
32 LauCartesianCPCoeffSet::LauCartesianCPCoeffSet(const TString& compName, Double_t x, Double_t y, Double_t deltaX, Double_t deltaY,
33  Bool_t xFixed, Bool_t yFixed, Bool_t deltaXFixed, Bool_t deltaYFixed, Bool_t deltaXSecondStage, Bool_t deltaYSecondStage) :
34  LauAbsCoeffSet(compName),
35  x_(new LauParameter("X", x, minRealImagPart_, maxRealImagPart_, xFixed)),
36  y_(new LauParameter("Y", y, minRealImagPart_, maxRealImagPart_, yFixed)),
37  deltaX_(new LauParameter("DeltaX", deltaX, minDelta_, maxDelta_, deltaXFixed)),
38  deltaY_(new LauParameter("DeltaY", deltaY, minDelta_, maxDelta_, deltaYFixed)),
39  particleCoeff_( x+deltaX, y+deltaY ),
40  antiparticleCoeff_( x-deltaX, y-deltaY ),
41  acp_("ACP", -2.0*(x*deltaX + y*deltaY)/(x*x + deltaX*deltaX + y*y + deltaY*deltaY), -1.0, 1.0, deltaXFixed&&deltaYFixed)
42 {
43  if (deltaXSecondStage && !deltaXFixed) {
44  deltaX_->secondStage(kTRUE);
45  deltaX_->initValue(0.0);
46  }
47  if (deltaYSecondStage && !deltaYFixed) {
48  deltaY_->secondStage(kTRUE);
49  deltaY_->initValue(0.0);
50  }
51 }
52 
54  x_(0),
55  y_(0),
56  deltaX_(0),
57  deltaY_(0),
58  particleCoeff_( rhs.particleCoeff_ ),
59  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
60  acp_( rhs.acp_ )
61 {
62  if ( cloneOption == All || cloneOption == TieRealPart ) {
63  x_ = rhs.x_->createClone(constFactor);
64  } else {
65  x_ = new LauParameter("X", rhs.x_->value(), minRealImagPart_, maxRealImagPart_, rhs.x_->fixed());
66  if ( rhs.x_->blind() ) {
67  const LauBlind* blinder = rhs.x_->blinder();
68  x_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
69  }
70  }
71 
72  if ( cloneOption == All || cloneOption == TieImagPart ) {
73  y_ = rhs.y_->createClone(constFactor);
74  } else {
75  y_ = new LauParameter("Y", rhs.y_->value(), minRealImagPart_, maxRealImagPart_, rhs.y_->fixed());
76  if ( rhs.y_->blind() ) {
77  const LauBlind* blinder = rhs.y_->blinder();
78  y_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
79  }
80  }
81 
82  if ( cloneOption == All || cloneOption == TieCPPars ) {
83  deltaX_ = rhs.deltaX_->createClone(constFactor);
84  deltaY_ = rhs.deltaY_->createClone(constFactor);
85  } else {
86  deltaX_ = new LauParameter("DeltaX", rhs.deltaX_->value(), minDelta_, maxDelta_, rhs.deltaX_->fixed());
87  deltaY_ = new LauParameter("DeltaY", rhs.deltaY_->value(), minDelta_, maxDelta_, rhs.deltaY_->fixed());
88  if ( rhs.deltaX_->secondStage() && !rhs.deltaX_->fixed() ) {
89  deltaX_->secondStage(kTRUE);
90  deltaX_->initValue(0.0);
91  }
92  if ( rhs.deltaY_->secondStage() && !rhs.deltaY_->fixed() ) {
93  deltaY_->secondStage(kTRUE);
94  deltaY_->initValue(0.0);
95  }
96  if ( rhs.deltaX_->blind() ) {
97  const LauBlind* blinder = rhs.deltaX_->blinder();
98  deltaX_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
99  }
100  if ( rhs.deltaY_->blind() ) {
101  const LauBlind* blinder = rhs.deltaY_->blinder();
102  deltaY_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
103  }
104  }
105 }
106 
107 std::vector<LauParameter*> LauCartesianCPCoeffSet::getParameters()
108 {
109  std::vector<LauParameter*> pars;
110  pars.push_back(x_);
111  pars.push_back(y_);
112  pars.push_back(deltaX_);
113  pars.push_back(deltaY_);
114  return pars;
115 }
116 
118 {
119  std::cout<<"INFO in LauCartesianCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
120  std::cout<<"x = "<<x_->value()<<",\t";
121  std::cout<<"y = "<<y_->value()<<",\t";
122  std::cout<<"Delta x = "<<deltaX_->value()<<",\t";
123  std::cout<<"Delta y = "<<deltaY_->value()<<"."<<std::endl;
124 }
125 
126 void LauCartesianCPCoeffSet::printTableHeading(std::ostream& stream) const
127 {
128  stream<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
129  stream<<"\\hline"<<std::endl;
130  stream<<"Component & Real Part & Imaginary Part & $\\Delta$ Real Part & $\\Delta$ Imaginary Part \\\\"<<std::endl;
131  stream<<"\\hline"<<std::endl;
132 }
133 
134 void LauCartesianCPCoeffSet::printTableRow(std::ostream& stream) const
135 {
136  LauPrint print;
137  TString resName = this->name();
138  resName = resName.ReplaceAll("_", "\\_");
139  stream<<resName<<" & $";
140  print.printFormat(stream, x_->value());
141  stream<<" \\pm ";
142  print.printFormat(stream, x_->error());
143  stream<<"$ & $";
144  print.printFormat(stream, y_->value());
145  stream<<" \\pm ";
146  print.printFormat(stream, y_->error());
147  stream<<"$ & $";
148  print.printFormat(stream, deltaX_->value());
149  stream<<" \\pm ";
150  print.printFormat(stream, deltaX_->error());
151  stream<<"$ & $";
152  print.printFormat(stream, deltaY_->value());
153  stream<<" \\pm ";
154  print.printFormat(stream, deltaY_->error());
155  stream<<"$ \\\\"<<std::endl;
156 }
157 
159 {
160  if (x_->fixed() == kFALSE) {
161  // Choose a value for "X" between -3.0 and 3.0
162  Double_t value = LauRandom::zeroSeedRandom()->Rndm()*6.0 - 3.0;
163  x_->initValue(value); x_->value(value);
164  }
165  if (y_->fixed() == kFALSE) {
166  // Choose a value for "Y" between -3.0 and 3.0
167  Double_t value = LauRandom::zeroSeedRandom()->Rndm()*6.0 - 3.0;
168  y_->initValue(value); y_->value(value);
169  }
170  if (deltaX_->fixed() == kFALSE && deltaX_->secondStage() == kFALSE) {
171  // Choose a value for "Delta X" between -0.5 and 0.5
172  Double_t value = LauRandom::zeroSeedRandom()->Rndm()*1.0 - 0.5;
173  deltaX_->initValue(value); deltaX_->value(value);
174  }
175  if (deltaY_->fixed() == kFALSE && deltaY_->secondStage() == kFALSE) {
176  // Choose a value for "Delta Y" between -0.5 and 0.5
177  Double_t value = LauRandom::zeroSeedRandom()->Rndm()*1.0 - 0.5;
178  deltaY_->initValue(value); deltaY_->value(value);
179  }
180 }
181 
183 {
184  // update the pulls
185  x_->updatePull();
186  y_->updatePull();
187  deltaX_->updatePull();
188  deltaY_->updatePull();
189 }
190 
192 {
194  return particleCoeff_;
195 }
196 
198 {
200  return antiparticleCoeff_;
201 }
202 
203 void LauCartesianCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
204 {
205  LauComplex average( coeff );
206  average += coeffBar;
207  average.rescale( 0.5 );
208 
209  Double_t xVal( average.re() );
210  Double_t yVal( average.im() );
211  Double_t deltaXVal( coeff.re() - average.re() );
212  Double_t deltaYVal( coeff.im() - average.im() );
213 
214  x_->value( xVal );
215  y_->value( yVal );
216  deltaX_->value( deltaXVal );
217  deltaY_->value( deltaYVal );
218 
219  if ( init ) {
220  x_->genValue( xVal );
221  y_->genValue( yVal );
222  deltaX_->genValue( deltaXVal );
223  deltaY_->genValue( deltaYVal );
224 
225  x_->initValue( xVal );
226  y_->initValue( yVal );
227  deltaX_->initValue( deltaXVal );
228  deltaY_->initValue( deltaYVal );
229  }
230 }
231 
233 {
234  // set the name
235  TString parName(this->baseName()); parName += "_ACP";
236  acp_.name(parName);
237 
238  // work out the ACP value
239  Double_t numer = x_->value()*deltaX_->value() + y_->value()*deltaY_->value();
240  Double_t denom = x_->value()*x_->value() + deltaX_->value()*deltaX_->value() + y_->value()*y_->value() + deltaY_->value()*deltaY_->value();
241  Double_t value = -2.0*numer/denom;
242 
243  // is it fixed?
244  Bool_t fixed = deltaX_->fixed() && deltaY_->fixed();
245  acp_.fixed(fixed);
246 
247  // we can't work out the error without the covariance matrix
248  Double_t error(0.0);
249 
250  // set the value and error
251  acp_.valueAndErrors(value,error);
252 
253  return acp_;
254 }
255 
256 LauAbsCoeffSet* LauCartesianCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
257 {
258  LauAbsCoeffSet* clone(0);
259  if ( cloneOption == All || cloneOption == TieRealPart || cloneOption == TieImagPart || cloneOption == TieCPPars ) {
260  clone = new LauCartesianCPCoeffSet( *this, cloneOption, constFactor );
261  clone->name( newName );
262  } else {
263  std::cerr << "ERROR in LauCartesianCPCoeffSet::createClone : Invalid clone option" << std::endl;
264  }
265  return clone;
266 }
267 
LauParameter * deltaY_
The asymmetric imaginary part.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
virtual void printParValues() const
Print the current values of the parameters.
static Double_t minDelta_
Minimum allowed value of CP-violating real/imaginary part parameters.
TRandom * zeroSeedRandom()
Access the singleton random number generator with seed set from machine clock time (within +-1 sec)...
Definition: LauRandom.cc:30
LauCartesianCPCoeffSet(const TString &compName, Double_t x, Double_t y, Double_t deltaX, Double_t deltaY, Bool_t xFixed, Bool_t yFixed, Bool_t deltaXFixed, Bool_t deltaYFixed, Bool_t deltaXSecondStage=kFALSE, Bool_t deltaYSecondStage=kFALSE)
Constructor.
ClassImp(LauAbsCoeffSet)
LauParameter * deltaX_
The asymmetric real part.
virtual void setCoeffValues(const LauComplex &coeff, const LauComplex &coeffBar, Bool_t init)
Set the parameters based on the complex coefficients for particles and antiparticles.
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
const TString & name() const
The parameter name.
LauParameter * x_
The average real part.
Double_t re() const
Get the real part.
Definition: LauComplex.hh:205
virtual LauParameter acp()
Calculate the CP asymmetry.
Double_t im() const
Get the imaginary part.
Definition: LauComplex.hh:214
static Double_t maxRealImagPart_
Maximum allowed value of real/imaginary part parameters.
LauComplex particleCoeff_
The particle complex coefficient.
static Double_t maxDelta_
Maximum allowed value of CP-violating real/imaginary part parameters.
File containing declaration of LauPrint class.
Class to define various output print commands.
Definition: LauPrint.hh:29
virtual const LauComplex & antiparticleCoeff()
Retrieve the complex coefficient for an antiparticle.
CloneOption
Options for cloning operation.
Bool_t clone() const
Check whether is a clone or not.
Bool_t blind() const
The blinding state.
virtual LauAbsCoeffSet * createClone(const TString &newName, CloneOption cloneOption=All, Double_t constFactor=1.0)
Create a clone of the coefficient set.
virtual void printTableRow(std::ostream &stream) const
Print the parameters of the complex coefficient as a row in the results table.
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
File containing declaration of LauParameter class.
LauParameter * y_
The average imaginary part.
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.
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.
static Double_t minRealImagPart_
Minimum allowed value of real/imaginary part parameters.
File containing declaration of LauCartesianCPCoeffSet class.
File containing LauRandom namespace.
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
void rescale(Double_t scaleVal)
Scale this by a factor.
Definition: LauComplex.hh:285
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.
Class for defining a complex number.
Definition: LauComplex.hh:47
void updatePull()
Call to update the bias and pull values.
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.
Double_t value() const
The value of 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 acp_
The CP asymmetry.
virtual void printTableHeading(std::ostream &stream) const
Print the column headings for a results table.
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.
Class for blinding and unblinding a number based on a blinding string.
Definition: LauBlind.hh:28
Double_t genValue() const
The value generated for the parameter.
Class for defining a complex coefficient using the Cartesian CP convention.
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.
virtual void randomiseInitValues()
Randomise the starting values of the parameters for a fit.