laura is hosted by Hepforge, IPPP Durham
Laura++  v3r4
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauCartesianGammaCPCoeffSet.cc
Go to the documentation of this file.
1 
2 /*
3 Copyright 2014 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 
37 #include "LauComplex.hh"
38 #include "LauConstants.hh"
39 #include "LauParameter.hh"
40 #include "LauPrint.hh"
41 
43 
44 
45 LauCartesianGammaCPCoeffSet::LauCartesianGammaCPCoeffSet(const TString& compName, const Double_t x, const Double_t y, const Double_t xCP, const Double_t yCP, const Double_t deltaXCP, const Double_t deltaYCP,
46  const Bool_t xFixed, const Bool_t yFixed, const Bool_t xCPFixed, const Bool_t yCPFixed, const Bool_t deltaXCPFixed, const Bool_t deltaYCPFixed, const Bool_t deltaXCPSecondStage, const Bool_t deltaYCPSecondStage) :
47  LauAbsCoeffSet(compName),
48  x_(new LauParameter("X", x, minRealImagPart_, maxRealImagPart_, xFixed)),
49  y_(new LauParameter("Y", y, minRealImagPart_, maxRealImagPart_, yFixed)),
50  xCP_(new LauParameter("XCP", xCP, minRealImagPart_, maxRealImagPart_, xCPFixed)),
51  yCP_(new LauParameter("YCP", yCP, minRealImagPart_, maxRealImagPart_, yCPFixed)),
52  deltaXCP_(new LauParameter("DeltaXCP", deltaXCP, minDelta_, maxDelta_, deltaXCPFixed)),
53  deltaYCP_(new LauParameter("DeltaYCP", deltaYCP, minDelta_, maxDelta_, deltaYCPFixed)),
54  nonCPPart_( x, y),
55  cpPart_( 1+xCP+deltaXCP, yCP+deltaYCP),
56  cpAntiPart_( 1+xCP-deltaXCP, yCP-deltaYCP),
57  particleCoeff_( nonCPPart_ * cpPart_ ),
58  antiparticleCoeff_( nonCPPart_ * cpAntiPart_ ),
59  acp_("ACP", (antiparticleCoeff_.abs2()-particleCoeff_.abs2())/(antiparticleCoeff_.abs2()+particleCoeff_.abs2()), -1.0, 1.0, deltaXCPFixed&&deltaYCPFixed)
60 {
61  if (deltaXCPSecondStage && !deltaXCPFixed) {
62  deltaXCP_->secondStage(kTRUE);
63  deltaXCP_->initValue(0.0);
64  }
65  if (deltaYCPSecondStage && !deltaYCPFixed) {
66  deltaYCP_->secondStage(kTRUE);
67  deltaYCP_->initValue(0.0);
68  }
69 }
70 
72  x_(0),
73  y_(0),
74  xCP_(0),
75  yCP_(0),
76  deltaXCP_(0),
77  deltaYCP_(0),
78  nonCPPart_( rhs.nonCPPart_ ),
79  cpPart_( rhs.cpPart_ ),
80  cpAntiPart_( rhs.cpAntiPart_ ),
81  particleCoeff_( rhs.particleCoeff_ ),
82  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
83  acp_( rhs.acp_ )
84 {
85  if ( cloneOption == All || cloneOption == TieRealPart ) {
86  x_ = rhs.x_->createClone(constFactor);
87  } else {
88  x_ = new LauParameter("X", rhs.x_->value(), minRealImagPart_, maxRealImagPart_, rhs.x_->fixed());
89  if ( rhs.x_->blind() ) {
90  const LauBlind* blinder = rhs.x_->blinder();
91  x_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
92  }
93  }
94 
95  if ( cloneOption == All || cloneOption == TieImagPart ) {
96  y_ = rhs.y_->createClone(constFactor);
97  } else {
98  y_ = new LauParameter("Y", rhs.y_->value(), minRealImagPart_, maxRealImagPart_, rhs.y_->fixed());
99  if ( rhs.y_->blind() ) {
100  const LauBlind* blinder = rhs.y_->blinder();
101  y_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
102  }
103  }
104 
105  if ( cloneOption == All || cloneOption == TieCPPars ) {
106  xCP_ = rhs.xCP_->createClone(constFactor);
107  yCP_ = rhs.yCP_->createClone(constFactor);
108  deltaXCP_ = rhs.deltaXCP_->createClone(constFactor);
109  deltaYCP_ = rhs.deltaYCP_->createClone(constFactor);
110  } else {
111  xCP_ = new LauParameter("XCP", rhs.xCP_->value(), minRealImagPart_, maxRealImagPart_, rhs.xCP_->fixed());
112  if ( rhs.xCP_->blind() ) {
113  const LauBlind* blinder = rhs.xCP_->blinder();
114  xCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
115  }
116  yCP_ = new LauParameter("YCP", rhs.yCP_->value(), minRealImagPart_, maxRealImagPart_, rhs.yCP_->fixed());
117  if ( rhs.yCP_->blind() ) {
118  const LauBlind* blinder = rhs.yCP_->blinder();
119  yCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
120  }
121  deltaXCP_ = new LauParameter("DeltaXCP", rhs.deltaXCP_->value(), minDelta_, maxDelta_, rhs.deltaXCP_->fixed());
122  deltaYCP_ = new LauParameter("DeltaYCP", rhs.deltaYCP_->value(), minDelta_, maxDelta_, rhs.deltaYCP_->fixed());
123  if ( rhs.deltaXCP_->secondStage() && !rhs.deltaXCP_->fixed() ) {
124  deltaXCP_->secondStage(kTRUE);
125  deltaXCP_->initValue(0.0);
126  }
127  if ( rhs.deltaYCP_->secondStage() && !rhs.deltaYCP_->fixed() ) {
128  deltaYCP_->secondStage(kTRUE);
129  deltaYCP_->initValue(0.0);
130  }
131  if ( rhs.deltaXCP_->blind() ) {
132  const LauBlind* blinder = rhs.deltaXCP_->blinder();
133  deltaXCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
134  }
135  if ( rhs.deltaYCP_->blind() ) {
136  const LauBlind* blinder = rhs.deltaYCP_->blinder();
137  deltaYCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
138  }
139  }
140 }
141 
142 std::vector<LauParameter*> LauCartesianGammaCPCoeffSet::getParameters()
143 {
144  std::vector<LauParameter*> pars;
145  pars.push_back(x_);
146  pars.push_back(y_);
147  if(!xCP_->fixed()) pars.push_back(xCP_);
148  if(!yCP_->fixed()) pars.push_back(yCP_);
149  if(!deltaXCP_->fixed()) pars.push_back(deltaXCP_);
150  if(!deltaYCP_->fixed()) pars.push_back(deltaYCP_);
151  return pars;
152 }
153 
155 {
156  std::cout<<"INFO in LauCartesianGammaCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
157  std::cout<<"x = "<<x_->value()<<",\t";
158  std::cout<<"y = "<<y_->value()<<",\t";
159  std::cout<<"xCP = "<<xCP_->value()<<",\t";
160  std::cout<<"yCP = "<<yCP_->value()<<",\t";
161  std::cout<<"Delta xCP = "<<deltaXCP_->value()<<",\t";
162  std::cout<<"Delta yCP = "<<deltaYCP_->value()<<"."<<std::endl;
163 }
164 
165 void LauCartesianGammaCPCoeffSet::printTableHeading(std::ostream& stream) const
166 {
167  stream<<"\\begin{tabular}{|l|c|c|c|c|c|c|}"<<std::endl;
168  stream<<"\\hline"<<std::endl;
169  stream<<"Component & Real Part & Imaginary Part & Real CP Part & Imaginary CP Part & $\\Delta$ Real CP Part & $\\Delta$ Imaginary CP Part \\\\"<<std::endl;
170  stream<<"\\hline"<<std::endl;
171 }
172 
173 void LauCartesianGammaCPCoeffSet::printTableRow(std::ostream& stream) const
174 {
175  LauPrint print;
176  TString resName = this->name();
177  resName = resName.ReplaceAll("_", "\\_");
178  stream<<resName<<" & $";
179  print.printFormat(stream, x_->value());
180  stream<<" \\pm ";
181  print.printFormat(stream, x_->error());
182  stream<<"$ & $";
183  print.printFormat(stream, y_->value());
184  stream<<" \\pm ";
185  print.printFormat(stream, y_->error());
186  stream<<"$ & $";
187  print.printFormat(stream, xCP_->value());
188  stream<<" \\pm ";
189  print.printFormat(stream, xCP_->error());
190  stream<<"$ & $";
191  print.printFormat(stream, yCP_->value());
192  stream<<" \\pm ";
193  print.printFormat(stream, yCP_->error());
194  stream<<"$ & $";
195  print.printFormat(stream, deltaXCP_->value());
196  stream<<" \\pm ";
197  print.printFormat(stream, deltaXCP_->error());
198  stream<<"$ & $";
199  print.printFormat(stream, deltaYCP_->value());
200  stream<<" \\pm ";
201  print.printFormat(stream, deltaYCP_->error());
202  stream<<"$ \\\\"<<std::endl;
203 }
204 
206 {
207  if (x_->fixed() == kFALSE) {
208  // Choose a value for "X" between -3.0 and 3.0
209  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
210  x_->initValue(value); x_->value(value);
211  }
212  if (y_->fixed() == kFALSE) {
213  // Choose a value for "Y" between -3.0 and 3.0
214  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
215  y_->initValue(value); y_->value(value);
216  }
217  if (xCP_->fixed() == kFALSE) {
218  // Choose a value for "XCP" between -3.0 and 3.0
219  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
220  xCP_->initValue(value); xCP_->value(value);
221  }
222  if (yCP_->fixed() == kFALSE) {
223  // Choose a value for "YCP" between -3.0 and 3.0
224  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
225  yCP_->initValue(value); yCP_->value(value);
226  }
227  if (deltaXCP_->fixed() == kFALSE && deltaXCP_->secondStage() == kFALSE) {
228  // Choose a value for "Delta XCP" between -0.5 and 0.5
229  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*1.0 - 0.5;
230  deltaXCP_->initValue(value); deltaXCP_->value(value);
231  }
232  if (deltaYCP_->fixed() == kFALSE && deltaYCP_->secondStage() == kFALSE) {
233  // Choose a value for "Delta YCP" between -0.5 and 0.5
234  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*1.0 - 0.5;
235  deltaYCP_->initValue(value); deltaYCP_->value(value);
236  }
237 }
238 
240 {
241  // update the pulls
242  x_->updatePull();
243  y_->updatePull();
244  xCP_->updatePull();
245  yCP_->updatePull();
248 }
249 
251 {
255  return particleCoeff_;
256 }
257 
259 {
263  return antiparticleCoeff_;
264 }
265 
267 {
268  std::cerr << "ERROR in LauCartesianGammaCPCoeffSet::setCoeffValues : Method not supported by this class - too many parameters" << std::endl;
269 }
270 
272 {
273  // set the name
274  TString parName(this->baseName()); parName += "_ACP";
275  acp_.name(parName);
276 
277  // work out the ACP value
278  const LauComplex nonCPPart( x_->value(), y_->value() );
279  const LauComplex cpPart( 1.0+xCP_->value()+deltaXCP_->value(), yCP_->value()+deltaYCP_->value() );
280  const LauComplex cpAntiPart( 1.0+xCP_->value()-deltaXCP_->value(), yCP_->value()-deltaYCP_->value() );
281  const LauComplex partCoeff = nonCPPart * cpPart;
282  const LauComplex antiCoeff = nonCPPart * cpAntiPart;
283 
284  const Double_t numer = antiCoeff.abs2() - partCoeff.abs2();
285  const Double_t denom = antiCoeff.abs2() + partCoeff.abs2();
286  const Double_t value = numer/denom;
287 
288  // is it fixed?
289  const Bool_t fixed = deltaXCP_->fixed() && deltaYCP_->fixed();
290  acp_.fixed(fixed);
291 
292  // we can't work out the error without the covariance matrix
293  const Double_t error(0.0);
294 
295  // set the value and error
296  acp_.valueAndErrors(value,error);
297 
298  return acp_;
299 }
300 
301 LauAbsCoeffSet* LauCartesianGammaCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
302 {
303  LauAbsCoeffSet* clone(0);
304  if ( cloneOption == All || cloneOption == TieRealPart || cloneOption == TieImagPart || cloneOption == TieCPPars ) {
305  clone = new LauCartesianGammaCPCoeffSet( *this, cloneOption, constFactor );
306  clone->name( newName );
307  } else {
308  std::cerr << "ERROR in LauCartesianGammaCPCoeffSet::createClone : Invalid clone option" << std::endl;
309  }
310  return clone;
311 }
312 
LauParameter * yCP_
The average CP imaginary part.
virtual LauParameter acp()
Calculate the CP asymmetry.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
LauComplex cpAntiPart_
The CP part of the complex coefficient for the antiparticle.
static Double_t minDelta_
Minimum allowed value of CP-violating real/imaginary part parameters.
virtual void printTableRow(std::ostream &stream) const
Print the parameters of the complex coefficient as a row in the results table.
virtual LauAbsCoeffSet * createClone(const TString &newName, CloneOption cloneOption=All, Double_t constFactor=1.0)
Create a clone of the coefficient set.
Class for defining a complex coefficient using the Cartesian gamma CP convention. ...
ClassImp(LauAbsCoeffSet)
LauParameter()
Default constructor.
Definition: LauParameter.cc:44
virtual void printTableHeading(std::ostream &stream) const
Print the column headings for a results table.
const TString & name() const
The parameter name.
LauComplex cpPart_
The CP part of the complex coefficient for the particle.
static Double_t maxRealImagPart_
Maximum allowed value of real/imaginary part parameters.
static Double_t maxDelta_
Maximum allowed value of CP-violating real/imaginary part parameters.
LauParameter * y_
The nonCP imaginary part.
File containing declaration of LauPrint class.
Class to define various output print commands.
Definition: LauPrint.hh:43
virtual std::vector< LauParameter * > getParameters()
Retrieve the parameters of the coefficient, e.g. so that they can be loaded into a fit...
CloneOption
Options for cloning operation.
Double_t abs2() const
Obtain the square of the absolute value of the complex number.
Definition: LauComplex.hh:246
Bool_t clone() const
Check whether is a clone or not.
virtual void randomiseInitValues()
Randomise the starting values of the parameters for a fit.
Bool_t blind() const
The blinding state.
LauComplex particleCoeff_
The particle complex coefficient.
virtual const LauComplex & antiparticleCoeff()
Retrieve the complex coefficient for an antiparticle.
LauParameter acp_
The CP asymmetry.
File containing declaration of LauParameter class.
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
Bool_t secondStage() const
Check whether the parameter should be floated only in the second stage of a two stage fit...
LauParameter * deltaXCP_
The asymmetric CP real part.
File containing declaration of LauComplex class.
LauParameter * x_
The nonCP real part.
const TString & blindingString() const
Obtain the blinding string.
Definition: LauBlind.hh:76
Double_t error() const
The error on the parameter.
LauParameter * xCP_
The average CP real part.
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 finaliseValues()
Make sure values are in &quot;standard&quot; ranges, e.g. phases should be between -pi and pi.
static Double_t minRealImagPart_
Minimum allowed value of real/imaginary part parameters.
const LauBlind * blinder() const
Access the blinder object.
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 * deltaYCP_
The asymmetric CP imaginary part.
void setRealImagPart(Double_t realpart, Double_t imagpart)
Set both real and imaginary part.
Definition: LauComplex.hh:328
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.
Class for defining a complex number.
Definition: LauComplex.hh:61
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.
File containing declaration of LauCartesianGammaCPCoeffSet class.
LauCartesianGammaCPCoeffSet(const TString &compName, const Double_t x, const Double_t y, const Double_t xCP, const Double_t yCP, const Double_t deltaXCP, const Double_t deltaYCP, const Bool_t xFixed, const Bool_t yFixed, const Bool_t xCPFixed, const Bool_t yCPFixed, const Bool_t deltaXCPFixed, const Bool_t deltaYCPFixed, const Bool_t deltaXCPSecondStage=kFALSE, const Bool_t deltaYCPSecondStage=kFALSE)
Constructor.
Double_t value() const
The value of the parameter.
virtual void printParValues() const
Print the current values of the parameters.
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.
LauComplex nonCPPart_
The nonCP part of the complex coefficient.
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.