laura is hosted by Hepforge, IPPP Durham
Laura++  v3r4
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauRealImagGammaCPCoeffSet.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 LauRealImagGammaCPCoeffSet::LauRealImagGammaCPCoeffSet(const TString& compName, const Double_t x, const Double_t y, const Double_t xCP, const Double_t yCP, const Double_t xbarCP, const Double_t ybarCP,
46  const Bool_t xFixed, const Bool_t yFixed, const Bool_t xCPFixed, const Bool_t yCPFixed, const Bool_t xbarCPFixed, const Bool_t ybarCPFixed) :
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  xbarCP_(new LauParameter("XbarCP", xbarCP, minRealImagPart_, maxRealImagPart_, xbarCPFixed)),
53  ybarCP_(new LauParameter("YbarCP", ybarCP, minRealImagPart_, maxRealImagPart_, ybarCPFixed)),
54  nonCPPart_( x, y),
55  cpPart_( 1+xCP, yCP ),
56  cpAntiPart_( 1+xbarCP, ybarCP ),
57  particleCoeff_( nonCPPart_ * cpPart_ ),
58  antiparticleCoeff_( nonCPPart_ * cpAntiPart_ ),
59  acp_("ACP", (antiparticleCoeff_.abs2()-particleCoeff_.abs2())/(antiparticleCoeff_.abs2()+particleCoeff_.abs2()), -1.0, 1.0, xCPFixed&&yCPFixed&&xbarCPFixed&&ybarCPFixed)
60 {
61 }
62 
64  x_(0),
65  y_(0),
66  xCP_(0),
67  yCP_(0),
68  xbarCP_(0),
69  ybarCP_(0),
70  nonCPPart_( rhs.nonCPPart_ ),
71  cpPart_( rhs.cpPart_ ),
72  cpAntiPart_( rhs.cpAntiPart_ ),
73  particleCoeff_( rhs.particleCoeff_ ),
74  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
75  acp_( rhs.acp_ )
76 {
77  if ( cloneOption == All || cloneOption == TieRealPart ) {
78  x_ = rhs.x_->createClone(constFactor);
79  } else {
80  x_ = new LauParameter("X", rhs.x_->value(), minRealImagPart_, maxRealImagPart_, rhs.x_->fixed());
81  if ( rhs.x_->blind() ) {
82  const LauBlind* blinder = rhs.x_->blinder();
83  x_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
84  }
85  }
86 
87  if ( cloneOption == All || cloneOption == TieImagPart ) {
88  y_ = rhs.y_->createClone(constFactor);
89  } else {
90  y_ = new LauParameter("Y", rhs.y_->value(), minRealImagPart_, maxRealImagPart_, rhs.y_->fixed());
91  if ( rhs.y_->blind() ) {
92  const LauBlind* blinder = rhs.y_->blinder();
93  y_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
94  }
95  }
96 
97  if ( cloneOption == All || cloneOption == TieCPPars ) {
98  xCP_ = rhs.xCP_->createClone(constFactor);
99  yCP_ = rhs.yCP_->createClone(constFactor);
100  xbarCP_ = rhs.xbarCP_->createClone(constFactor);
101  ybarCP_ = rhs.ybarCP_->createClone(constFactor);
102  } else {
103  xCP_ = new LauParameter("XCP", rhs.xCP_->value(), minRealImagPart_, maxRealImagPart_, rhs.xCP_->fixed());
104  if ( rhs.xCP_->blind() ) {
105  const LauBlind* blinder = rhs.xCP_->blinder();
106  xCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
107  }
108  yCP_ = new LauParameter("YCP", rhs.yCP_->value(), minRealImagPart_, maxRealImagPart_, rhs.yCP_->fixed());
109  if ( rhs.yCP_->blind() ) {
110  const LauBlind* blinder = rhs.yCP_->blinder();
111  yCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
112  }
114  if ( rhs.xbarCP_->blind() ) {
115  const LauBlind* blinder = rhs.xbarCP_->blinder();
116  xbarCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
117  }
119  if ( rhs.ybarCP_->blind() ) {
120  const LauBlind* blinder = rhs.ybarCP_->blinder();
121  ybarCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
122  }
123  }
124 }
125 
126 std::vector<LauParameter*> LauRealImagGammaCPCoeffSet::getParameters()
127 {
128  std::vector<LauParameter*> pars;
129  pars.push_back(x_);
130  pars.push_back(y_);
131  if(!xCP_->fixed()) pars.push_back(xCP_);
132  if(!yCP_->fixed()) pars.push_back(yCP_);
133  if(!xbarCP_->fixed()) pars.push_back(xbarCP_);
134  if(!ybarCP_->fixed()) pars.push_back(ybarCP_);
135  return pars;
136 }
137 
139 {
140  std::cout << "INFO in LauRealImagGammaCPCoeffSet::printParValues : Component \"" << this->name() << "\" has ";
141  std::cout << "x = " << x_->value() << ",\t";
142  std::cout << "y = " << y_->value() << ",\t";
143  std::cout << "xCP = " << xCP_->value() << ",\t";
144  std::cout << "yCP = " << yCP_->value() << ",\t";
145  std::cout << "xbarCP = " << xbarCP_->value() << ",\t";
146  std::cout << "ybarCP = " << ybarCP_->value() << "." << std::endl;
147 }
148 
149 void LauRealImagGammaCPCoeffSet::printTableHeading(std::ostream& stream) const
150 {
151  stream<<"\\begin{tabular}{|l|c|c|c|c|c|c|}"<<std::endl;
152  stream<<"\\hline"<<std::endl;
153  stream<<"Component & Real Part & Imaginary Part & Particle CP Real Part & Particle CP Imaginary Part & Antiparticle CP Real Part & Antiparticle CP Imaginary Part \\\\"<<std::endl;
154  stream<<"\\hline"<<std::endl;
155 }
156 
157 void LauRealImagGammaCPCoeffSet::printTableRow(std::ostream& stream) const
158 {
159  LauPrint print;
160  TString resName = this->name();
161  resName = resName.ReplaceAll("_", "\\_");
162  stream<<resName<<" & $";
163  print.printFormat(stream, x_->value());
164  stream<<" \\pm ";
165  print.printFormat(stream, x_->error());
166  stream<<"$ & $";
167  print.printFormat(stream, y_->value());
168  stream<<" \\pm ";
169  print.printFormat(stream, y_->error());
170  stream<<"$ & $";
171  print.printFormat(stream, xCP_->value());
172  stream<<" \\pm ";
173  print.printFormat(stream, xCP_->error());
174  stream<<"$ & $";
175  print.printFormat(stream, yCP_->value());
176  stream<<" \\pm ";
177  print.printFormat(stream, yCP_->error());
178  stream<<"$ & $";
179  print.printFormat(stream, xbarCP_->value());
180  stream<<" \\pm ";
181  print.printFormat(stream, xbarCP_->error());
182  stream<<"$ & $";
183  print.printFormat(stream, ybarCP_->value());
184  stream<<" \\pm ";
185  print.printFormat(stream, ybarCP_->error());
186  stream<<"$ \\\\"<<std::endl;
187 }
188 
190 {
191  if (x_->fixed() == kFALSE) {
192  // Choose a value for "X" between -3.0 and 3.0
193  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
194  x_->initValue(value); x_->value(value);
195  }
196  if (y_->fixed() == kFALSE) {
197  // Choose a value for "Y" between -3.0 and 3.0
198  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
199  y_->initValue(value); y_->value(value);
200  }
201  if (xCP_->fixed() == kFALSE) {
202  // Choose a value for "XCP" between -3.0 and 3.0
203  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
204  xCP_->initValue(value); xCP_->value(value);
205  }
206  if (yCP_->fixed() == kFALSE) {
207  // Choose a value for "YCP" between -3.0 and 3.0
208  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
209  yCP_->initValue(value); yCP_->value(value);
210  }
211  if (xbarCP_->fixed() == kFALSE) {
212  // Choose a value for "XbarCP" between -3.0 and 3.0
213  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
214  xbarCP_->initValue(value); xbarCP_->value(value);
215  }
216  if (ybarCP_->fixed() == kFALSE) {
217  // Choose a value for "YbarCP" between -3.0 and 3.0
218  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
219  ybarCP_->initValue(value); ybarCP_->value(value);
220  }
221 }
222 
224 {
225  // update the pulls
226  x_->updatePull();
227  y_->updatePull();
228  xCP_->updatePull();
229  yCP_->updatePull();
230  xbarCP_->updatePull();
231  ybarCP_->updatePull();
232 }
233 
235 {
239  return particleCoeff_;
240 }
241 
243 {
247  return antiparticleCoeff_;
248 }
249 
251 {
252  std::cerr << "ERROR in LauCartesianGammaCPCoeffSet::setCoeffValues : Method not supported by this class - too many parameters" << std::endl;
253 }
254 
256 {
257  // set the name
258  TString parName(this->baseName()); parName += "_ACP";
259  acp_.name(parName);
260 
261  // work out the ACP value
262  const LauComplex nonCPPart( x_->value(), y_->value() );
263  const LauComplex cpPart( 1.0+xCP_->value(), yCP_->value() );
264  const LauComplex cpAntiPart( 1.0+xbarCP_->value(), ybarCP_->value() );
265  const LauComplex partCoeff = nonCPPart * cpPart;
266  const LauComplex antiCoeff = nonCPPart * cpAntiPart;
267 
268  const Double_t numer = antiCoeff.abs2() - partCoeff.abs2();
269  const Double_t denom = antiCoeff.abs2() + partCoeff.abs2();
270  const Double_t value = numer/denom;
271 
272  // is it fixed?
273  const Bool_t fixed = xCP_->fixed() && yCP_->fixed() && xbarCP_->fixed() && ybarCP_->fixed();
274  acp_.fixed(fixed);
275 
276  // we can't work out the error without the covariance matrix
277  const Double_t error(0.0);
278 
279  // set the value and error
280  acp_.valueAndErrors(value,error);
281 
282  return acp_;
283 }
284 
285 LauAbsCoeffSet* LauRealImagGammaCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
286 {
287  LauAbsCoeffSet* clone(0);
288  if ( cloneOption == All || cloneOption == TieRealPart || cloneOption == TieImagPart || cloneOption == TieCPPars ) {
289  clone = new LauRealImagGammaCPCoeffSet( *this, cloneOption, constFactor );
290  clone->name( newName );
291  } else {
292  std::cerr << "ERROR in LauRealImagGammaCPCoeffSet::createClone : Invalid clone option" << std::endl;
293  }
294  return clone;
295 }
296 
LauRealImagGammaCPCoeffSet(const TString &compName, const Double_t x, const Double_t y, const Double_t xCP, const Double_t yCP, const Double_t xbarCP, const Double_t ybarCP, const Bool_t xFixed, const Bool_t yFixed, const Bool_t xCPFixed, const Bool_t yCPFixed, const Bool_t xbarCPFixed, const Bool_t ybarCPFixed)
Constructor.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
virtual LauAbsCoeffSet * createClone(const TString &newName, CloneOption cloneOption=All, Double_t constFactor=1.0)
Create a clone of the coefficient set.
LauParameter * xbarCP_
The real CP part for the antiparticle.
LauComplex cpAntiPart_
The CP part of the complex coefficient for the antiparticle.
ClassImp(LauAbsCoeffSet)
virtual void printParValues() const
Print the current values of the parameters.
LauParameter()
Default constructor.
Definition: LauParameter.cc:44
const TString & name() const
The parameter name.
Class for defining a complex coefficient using a Cartesian nonCP part multiplied by a simple Cartesia...
static Double_t maxRealImagPart_
Maximum allowed value of 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.
File containing declaration of LauPrint class.
Class to define various output print commands.
Definition: LauPrint.hh:43
LauParameter * ybarCP_
The imaginary CP part for the antiparticle.
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
LauComplex particleCoeff_
The particle complex coefficient.
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.
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.
LauComplex cpPart_
The CP part of the complex coefficient for the particle.
Bool_t blind() const
The blinding state.
LauParameter * yCP_
The imaginary CP part for the particle.
LauParameter * xCP_
The real CP part for the particle.
File containing declaration of LauParameter class.
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.
Class for defining the abstract interface for complex coefficient classes.
virtual void setCoeffValues(const LauComplex &coeff, const LauComplex &coeffBar, Bool_t init)
Set the parameters based on the complex coefficients for particles and antiparticles.
File containing declaration of LauRealImagGammaCPCoeffSet class.
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 const LauComplex & antiparticleCoeff()
Retrieve the complex coefficient for an antiparticle.
virtual std::vector< LauParameter * > getParameters()
Retrieve the parameters of the coefficient, e.g. so that they can be loaded into a fit...
static Double_t minRealImagPart_
Minimum allowed value of real/imaginary part parameters.
LauParameter * y_
The imaginary nonCP part.
const LauBlind * blinder() const
Access the blinder object.
LauParameter acp_
The CP asymmetry.
void setRealImagPart(Double_t realpart, Double_t imagpart)
Set both real and imaginary part.
Definition: LauComplex.hh:328
virtual void printTableHeading(std::ostream &stream) const
Print the column headings for a results table.
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.
LauParameter * x_
The real nonCP part.
Double_t value() const
The value of the parameter.
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 randomiseInitValues()
Randomise the starting values of the parameters for a fit.
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.
virtual void finaliseValues()
Make sure values are in &quot;standard&quot; ranges, e.g. phases should be between -pi and pi.
virtual LauParameter acp()
Calculate the CP asymmetry.