laura is hosted by Hepforge, IPPP Durham
Laura++  v3r3
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauMagPhaseCPCoeffSet.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2011 - 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 "LauMagPhaseCPCoeffSet.hh"
23 #include "LauComplex.hh"
24 #include "LauConstants.hh"
25 #include "LauParameter.hh"
26 #include "LauPrint.hh"
27 
29 
30 
31 LauMagPhaseCPCoeffSet::LauMagPhaseCPCoeffSet(const TString& compName, Double_t mag, Double_t phase, Double_t magBar, Double_t phaseBar,
32  Bool_t magFixed, Bool_t phaseFixed,Bool_t magBarFixed, Bool_t phaseBarFixed) :
33  LauAbsCoeffSet(compName),
34  mag_(new LauParameter("Mag", mag, minMagnitude_, maxMagnitude_, magFixed)),
35  phase_(new LauParameter("Phase", phase, minPhase_, maxPhase_, phaseFixed)),
36  magBar_(new LauParameter("MagBar", magBar, minMagnitude_, maxMagnitude_, magBarFixed)),
37  phaseBar_(new LauParameter("PhaseBar", phaseBar, minPhase_, maxPhase_, phaseBarFixed)),
38  particleCoeff_( mag*TMath::Cos(phase), mag*TMath::Sin(phase) ),
39  antiparticleCoeff_( magBar*TMath::Cos(phaseBar), magBar*TMath::Sin(phaseBar) ),
40  acp_("ACP", (magBar*magBar - mag*mag)/(magBar*magBar + mag*mag), -1.0, 1.0)
41 {
42 }
43 
45  mag_(0),
46  phase_(0),
47  magBar_(0),
48  phaseBar_(0),
49  particleCoeff_( rhs.particleCoeff_ ),
50  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
51  acp_( rhs.acp_ )
52 {
53  if ( cloneOption == All || cloneOption == TieMagnitude ) {
54  mag_ = rhs.mag_->createClone(constFactor);
55  magBar_ = rhs.magBar_->createClone(constFactor);
56  } else {
57  mag_ = new LauParameter("Mag", rhs.mag_->value(), minMagnitude_, maxMagnitude_, rhs.mag_->fixed());
58  if ( rhs.mag_->blind() ) {
59  const LauBlind* blinder = rhs.mag_->blinder();
60  mag_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
61  }
62  magBar_ = new LauParameter("MagBar", rhs.magBar_->value(), minMagnitude_, maxMagnitude_, rhs.magBar_->fixed());
63  if ( rhs.magBar_->blind() ) {
64  const LauBlind* blinder = rhs.magBar_->blinder();
65  magBar_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
66  }
67  }
68 
69  if ( cloneOption == All || cloneOption == TiePhase ) {
70  phase_ = rhs.phase_->createClone(constFactor);
71  phaseBar_ = rhs.phaseBar_->createClone(constFactor);
72  } else {
73  phase_ = new LauParameter("Phase", rhs.phase_->value(), minPhase_, maxPhase_, rhs.phase_->fixed());
74  if ( rhs.phase_->blind() ) {
75  const LauBlind* blinder = rhs.phase_->blinder();
76  phase_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
77  }
78  phaseBar_ = new LauParameter("PhaseBar", rhs.phaseBar_->value(), minPhase_, maxPhase_, rhs.phaseBar_->fixed());
79  if ( rhs.phaseBar_->blind() ) {
80  const LauBlind* blinder = rhs.phaseBar_->blinder();
81  phaseBar_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
82  }
83  }
84 }
85 
86 std::vector<LauParameter*> LauMagPhaseCPCoeffSet::getParameters()
87 {
88  std::vector<LauParameter*> pars;
89  pars.push_back(mag_);
90  pars.push_back(phase_);
91  pars.push_back(magBar_);
92  pars.push_back(phaseBar_);
93  return pars;
94 }
95 
97 {
98  std::cout<<"INFO in LauMagPhaseCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
99  std::cout<<"mag = "<<mag_->value()<<",\t";
100  std::cout<<"phase = "<<phase_->value()<<",\t";
101  std::cout<<"magBar = "<<magBar_->value()<<",\t";
102  std::cout<<"phaseBar = "<<phaseBar_->value()<<"."<<std::endl;
103 }
104 
105 void LauMagPhaseCPCoeffSet::printTableHeading(std::ostream& stream) const
106 {
107  stream<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
108  stream<<"\\hline"<<std::endl;
109  stream<<"Component & Magnitude & Phase & Magnitude_bar & Phase_bar \\\\"<<std::endl;
110  stream<<"\\hline"<<std::endl;
111 }
112 
113 void LauMagPhaseCPCoeffSet::printTableRow(std::ostream& stream) const
114 {
115  LauPrint print;
116  TString resName = this->name();
117  resName = resName.ReplaceAll("_", "\\_");
118  stream<<resName<<" & $";
119  print.printFormat(stream, mag_->value());
120  stream<<" \\pm ";
121  print.printFormat(stream, mag_->error());
122  stream<<"$ & $";
123  print.printFormat(stream, phase_->value());
124  stream<<" \\pm ";
125  print.printFormat(stream, phase_->error());
126  stream<<"$ & $";
127  print.printFormat(stream, magBar_->value());
128  stream<<" \\pm ";
129  print.printFormat(stream, magBar_->error());
130  stream<<"$ & $";
131  print.printFormat(stream, phaseBar_->value());
132  stream<<" \\pm ";
133  print.printFormat(stream, phaseBar_->error());
134  stream<<"$ \\\\"<<std::endl;
135 }
136 
138 {
139  if (mag_->fixed() == kFALSE) {
140  // Choose a value for "magnitude" between 0.0 and 2.0
141  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*2.0;
142  mag_->initValue(value); mag_->value(value);
143  }
144  if (phase_->fixed() == kFALSE) {
145  // Choose a phase between +- pi
147  phase_->initValue(value); phase_->value(value);
148  }
149  if (magBar_->fixed() == kFALSE) {
150  // Choose a value for "magnitude" between 0.0 and 2.0
151  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*2.0;
152  magBar_->initValue(value); magBar_->value(value);
153  }
154  if (phaseBar_->fixed() == kFALSE) {
155  // Choose a phase between +- pi
157  phaseBar_->initValue(value); phaseBar_->value(value);
158  }
159 }
160 
162 {
163  // retrieve the current values from the parameters
164  Double_t mVal= mag_->value();
165  Double_t pVal= phase_->value();
166  Double_t mBarVal= magBar_->value();
167  Double_t pBarVal= phaseBar_->value();
168  Double_t genPhase = phase_->genValue();
169  Double_t genPhaseBar = phaseBar_->genValue();
170 
171  // Check whether we have a negative magnitude.
172  // If so make it positive and add pi to the phase.
173  if (mVal < 0.0) {
174  mVal *= -1.0;
175  pVal += LauConstants::pi;
176  }
177  if (mBarVal < 0.0) {
178  mBarVal *= -1.0;
179  pBarVal += LauConstants::pi;
180  }
181 
182  // Check now whether the phases lies in the right range (-pi to pi).
183  Bool_t pWithinRange(kFALSE);
184  Bool_t pBarWithinRange(kFALSE);
185  while (pWithinRange == kFALSE && pBarWithinRange == kFALSE) {
186  if (pVal > -LauConstants::pi && pVal < LauConstants::pi) {
187  pWithinRange = kTRUE;
188  } else {
189  // Not within the specified range
190  if (pVal > LauConstants::pi) {
191  pVal -= LauConstants::twoPi;
192  } else if (pVal < -LauConstants::pi) {
193  pVal += LauConstants::twoPi;
194  }
195  }
196 
197  if (pBarVal > -LauConstants::pi && pBarVal < LauConstants::pi) {
198  pBarWithinRange = kTRUE;
199  } else {
200  // Not within the specified range
201  if (pBarVal > LauConstants::pi) {
202  pBarVal -= LauConstants::twoPi;
203  } else if (pBarVal < -LauConstants::pi) {
204  pBarVal += LauConstants::twoPi;
205  }
206  }
207  }
208 
209  // A further problem can occur when the generated phase is close to -pi or pi.
210  // The phase can wrap over to the other end of the scale -
211  // this leads to artificially large pulls so we wrap it back.
212  Double_t diff = pVal - genPhase;
213  if (diff > LauConstants::pi) {
214  pVal -= LauConstants::twoPi;
215  } else if (diff < -LauConstants::pi) {
216  pVal += LauConstants::twoPi;
217  }
218 
219  diff = pBarVal - genPhaseBar;
220  if (diff > LauConstants::pi) {
221  pBarVal -= LauConstants::twoPi;
222  } else if (diff < -LauConstants::pi) {
223  pBarVal += LauConstants::twoPi;
224  }
225 
226  // finally store the new values in the parameters
227  // and update the pulls
228  mag_->value(mVal); mag_->updatePull();
229  phase_->value(pVal); phase_->updatePull();
230  magBar_->value(mBarVal); magBar_->updatePull();
231  phaseBar_->value(pBarVal); phaseBar_->updatePull();
232 
233 
234 }
235 
237 {
239  return particleCoeff_;
240 }
241 
243 {
245  return antiparticleCoeff_;
246 }
247 
248 void LauMagPhaseCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
249 {
250  Double_t magVal( coeff.abs() );
251  Double_t phaseVal( coeff.arg() );
252  Double_t magBarVal( coeffBar.abs() );
253  Double_t phaseBarVal( coeffBar.arg() );
254 
255  mag_->value( magVal );
256  phase_->value( phaseVal );
257  magBar_->value( magBarVal );
258  phaseBar_->value( phaseBarVal );
259 
260  if ( init ) {
261  mag_->genValue( magVal );
262  phase_->genValue( phaseVal );
263  magBar_->genValue( magBarVal );
264  phaseBar_->genValue( phaseBarVal );
265 
266  mag_->initValue( magVal );
267  phase_->initValue( phaseVal );
268  magBar_->initValue( magBarVal );
269  phaseBar_->initValue( phaseBarVal );
270  }
271 }
272 
274 {
275  // set the name
276  TString parName(this->baseName()); parName += "_ACP";
277  acp_.name(parName);
278 
279  // work out the ACP value
280 
281  Double_t value(-99.0);
282  value = (magBar_->value()*magBar_->value() - mag_->value()*mag_->value())/(magBar_->value()*magBar_->value() + mag_->value()*mag_->value());
283 
284  // is it fixed?
285  Bool_t fixed = magBar_->fixed() && mag_->fixed();
286  acp_.fixed(fixed);
287 
288  // we can't work out the error without the covariance matrix
289  Double_t error(0.0);
290 
291  // set the value and error
292  acp_.valueAndErrors(value,error);
293 
294  return acp_;
295 }
296 
297 LauAbsCoeffSet* LauMagPhaseCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
298 {
299  LauAbsCoeffSet* clone(0);
300  if ( cloneOption == All || cloneOption == TiePhase || cloneOption == TieMagnitude ) {
301  clone = new LauMagPhaseCPCoeffSet( *this, cloneOption, constFactor );
302  clone->name( newName );
303  } else {
304  std::cerr << "ERROR in LauMagPhaseCPCoeffSet::createClone : Invalid clone option" << std::endl;
305  }
306  return clone;
307 }
308 
static Double_t maxPhase_
Maximum allowed value of phase parameters.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
virtual void finaliseValues()
Make sure values are in &quot;standard&quot; ranges, e.g. phases should be between -pi and pi.
LauParameter acp_
The CP asymmetry.
virtual void randomiseInitValues()
Randomise the starting values of the parameters for a fit.
const Double_t twoPi
Two times Pi.
Definition: LauConstants.hh:93
File containing declaration of LauMagPhaseCPCoeffSet class.
ClassImp(LauAbsCoeffSet)
LauParameter * phase_
The phase for particles.
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
const TString & name() const
The parameter name.
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.
virtual LauAbsCoeffSet * createClone(const TString &newName, CloneOption cloneOption=All, Double_t constFactor=1.0)
Create a clone of the coefficient set.
virtual std::vector< LauParameter * > getParameters()
Retrieve the parameters of the coefficient, e.g. so that they can be loaded into a fit...
Class to define various output print commands.
Definition: LauPrint.hh:29
Class for defining a complex coefficient using seperate magnitudes and phases for particles and antip...
virtual void printParValues() const
Print the current values of the parameters.
CloneOption
Options for cloning operation.
Bool_t clone() const
Check whether is a clone or not.
Bool_t blind() const
The blinding state.
static Double_t maxMagnitude_
Maximum allowed value of magnitude parameters.
File containing declaration of LauParameter class.
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
File containing declaration of LauComplex class.
const TString & blindingString() const
Obtain the blinding string.
Definition: LauBlind.hh:62
virtual void setCoeffValues(const LauComplex &coeff, const LauComplex &coeffBar, Bool_t init)
Set the parameters based on the complex coefficients for particles and antiparticles.
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.
virtual LauParameter acp()
Calculate the CP asymmetry.
LauParameter * magBar_
The magnitude for antiparticles.
const LauBlind * blinder() const
Access the blinder object.
virtual const LauComplex & antiparticleCoeff()
Retrieve the complex coefficient for an antiparticle.
void setRealImagPart(Double_t realpart, Double_t imagpart)
Set both real and imaginary part.
Definition: LauComplex.hh:314
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 * mag_
The magnitude for particles.
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.
LauMagPhaseCPCoeffSet(const TString &compName, Double_t mag, Double_t phase, Double_t magBar, Double_t phaseBar, Bool_t magFixed, Bool_t phaseFixed, Bool_t magBarFixed, Bool_t phaseBarFixed)
Constructor.
static Double_t minPhase_
Minimum allowed value of phase parameters.
Double_t value() const
The value of the parameter.
LauComplex particleCoeff_
The particle complex coefficient.
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
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.
LauParameter * phaseBar_
The phase for antiparticles.
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.
virtual void printTableHeading(std::ostream &stream) const
Print the column headings for a results table.
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.