laura is hosted by Hepforge, IPPP Durham
Laura++  v3r1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauPolarGammaCPCoeffSet.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 
32 
34 
35 
36 LauPolarGammaCPCoeffSet::LauPolarGammaCPCoeffSet(const TString& compName, const DecayType decayType,
37  const Double_t x, const Double_t y,
38  const Double_t rB, const Double_t deltaB, const Double_t gamma,
39  const Double_t rD, const Double_t deltaD,
40  const Bool_t xFixed, const Bool_t yFixed,
41  const Bool_t rBFixed, const Bool_t deltaBFixed, const Bool_t gammaFixed,
42  const Bool_t rDFixed, const Bool_t deltaDFixed,
43  const Bool_t rBSecondStage, const Bool_t deltaBSecondStage, const Bool_t gammaSecondStage,
44  const Bool_t rDSecondStage, const Bool_t deltaDSecondStage,
45  const Bool_t useGlobalGamma,
46  const Bool_t useGlobalADSPars) :
47  LauAbsCoeffSet(compName),
48  decayType_(decayType),
49  x_(0),
50  y_(0),
51  rB_(0),
52  deltaB_(0),
53  gamma_(0),
54  rD_(0),
55  deltaD_(0),
56  useGlobalGamma_(useGlobalGamma),
57  useGlobalADSPars_(useGlobalADSPars),
58  nonCPPart_(x,y),
59  cpPart_(0.0,0.0),
60  cpAntiPart_(0.0,0.0),
61  particleCoeff_(0.0,0.0),
62  antiparticleCoeff_(0.0,0.0),
63  acp_("ACP", 0.0, -1.0, 1.0)
64 {
65  // All of the possible D decay types need these two parameters
66  x_ = new LauParameter("X", x, minRealImagPart_, maxRealImagPart_, xFixed);
67  y_ = new LauParameter("Y", y, minRealImagPart_, maxRealImagPart_, yFixed);
68 
69  // if we're using a global gamma, create it if it doesn't already exist then set gamma_ to point to it
70  // otherwise create our individual copy of gamma
71  if (useGlobalGamma_) {
72  if (!gammaGlobal_) {
73  gammaGlobal_ = new LauParameter("gamma", gamma, minPhase_, maxPhase_, gammaFixed);
74  gamma_ = gammaGlobal_;
75  } else {
76  gamma_ = gammaGlobal_->createClone();
77  }
78  } else {
79  gamma_ = new LauParameter("gamma", gamma, minPhase_, maxPhase_, gammaFixed);
80  }
81  if (gammaSecondStage && !gammaFixed) {
82  gamma_->secondStage(kTRUE);
83  gamma_->initValue(0.0);
84  }
85 
86  // which of the other parameter we need depends on the D-decay type
87  if ( decayType_ == ADS_Favoured || decayType_ == ADS_Suppressed || decayType_ == GLW_CPOdd || decayType_ == GLW_CPEven ) {
88  rB_ = new LauParameter("rB", rB, minMagnitude_, maxMagnitude_, rBFixed);
89  deltaB_ = new LauParameter("deltaB", deltaB, minPhase_, maxPhase_, deltaBFixed);
90  if (rBSecondStage && !rBFixed) {
91  rB_->secondStage(kTRUE);
92  rB_->initValue(0.0);
93  }
94  if (deltaBSecondStage && !deltaBFixed) {
95  deltaB_->secondStage(kTRUE);
96  deltaB_->initValue(0.0);
97  }
98  }
99 
100  if ( decayType_ == ADS_Favoured || decayType_ == ADS_Suppressed || decayType_ == ADS_Favoured_btouOnly ) {
101  if (useGlobalADSPars_) {
102  if ( !rDGlobal_ ) {
103  rDGlobal_ = new LauParameter("rD", rD, minMagnitude_, maxMagnitude_, rDFixed);
104  deltaDGlobal_ = new LauParameter("deltaD", deltaD, minPhase_, maxPhase_, deltaDFixed);
105  rD_ = rDGlobal_;
106  deltaD_ = deltaDGlobal_;
107  } else {
108  rD_ = rDGlobal_->createClone();
109  deltaD_ = deltaDGlobal_->createClone();
110  }
111  } else {
112  rD_ = new LauParameter("rD", rD, minMagnitude_, maxMagnitude_, rDFixed);
113  deltaD_ = new LauParameter("deltaD", deltaD, minPhase_, maxPhase_, deltaDFixed);
114  }
115  if (rDSecondStage && !rDFixed) {
116  rD_->secondStage(kTRUE);
117  rD_->initValue(0.0);
118  }
119  if (deltaDSecondStage && !deltaDFixed) {
120  deltaD_->secondStage(kTRUE);
121  deltaD_->initValue(0.0);
122  }
123  }
124 }
125 
127  decayType_( rhs.decayType_ ),
128  x_(0),
129  y_(0),
130  rB_(0),
131  deltaB_(0),
132  gamma_(0),
133  rD_(0),
134  deltaD_(0),
135  useGlobalGamma_( rhs.useGlobalGamma_ ),
136  useGlobalADSPars_( rhs.useGlobalADSPars_ ),
137  nonCPPart_( rhs.nonCPPart_ ),
138  cpPart_( rhs.cpPart_ ),
139  cpAntiPart_( rhs.cpAntiPart_ ),
140  particleCoeff_( rhs.particleCoeff_ ),
141  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
142  acp_( rhs.acp_ )
143 {
144  if ( cloneOption == All || cloneOption == TieRealPart ) {
145  x_ = rhs.x_->createClone(constFactor);
146  } else {
147  x_ = new LauParameter("X", rhs.x_->value(), minRealImagPart_, maxRealImagPart_, rhs.x_->fixed());
148  if ( rhs.x_->blind() ) {
149  const LauBlind* blinder = rhs.x_->blinder();
150  x_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
151  }
152  }
153 
154  if ( cloneOption == All || cloneOption == TieImagPart ) {
155  y_ = rhs.y_->createClone(constFactor);
156  } else {
157  y_ = new LauParameter("Y", rhs.y_->value(), minRealImagPart_, maxRealImagPart_, rhs.y_->fixed());
158  if ( rhs.y_->blind() ) {
159  const LauBlind* blinder = rhs.y_->blinder();
160  y_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
161  }
162  }
163 
164  if ( cloneOption == All || cloneOption == TieCPPars ) {
165  gamma_ = rhs.gamma_->createClone(constFactor);
167  rB_ = rhs.rB_->createClone(constFactor);
168  deltaB_ = rhs.deltaB_->createClone(constFactor);
169  }
171  rD_ = rhs.rD_->createClone(constFactor);
172  deltaD_ = rhs.deltaD_->createClone(constFactor);
173  }
174  } else {
175  if (useGlobalGamma_) {
177  } else {
178  gamma_ = new LauParameter("gamma", rhs.gamma_->value(), minPhase_, maxPhase_, rhs.gamma_->fixed());
179  if ( rhs.gamma_->blind() ) {
180  const LauBlind* blinder = rhs.gamma_->blinder();
181  gamma_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
182  }
183  if ( rhs.gamma_->secondStage() && !rhs.gamma_->fixed() ) {
184  gamma_->secondStage(kTRUE);
185  gamma_->initValue(0.0);
186  }
187  }
189  rB_ = new LauParameter("rB", rhs.rB_->value(), minMagnitude_, maxMagnitude_, rhs.rB_->fixed());
190  if ( rhs.rB_->blind() ) {
191  const LauBlind* blinder = rhs.rB_->blinder();
192  rB_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
193  }
194  deltaB_ = new LauParameter("deltaB", rhs.deltaB_->value(), minPhase_, maxPhase_, rhs.deltaB_->fixed());
195  if ( rhs.deltaB_->blind() ) {
196  const LauBlind* blinder = rhs.deltaB_->blinder();
197  deltaB_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
198  }
199  if ( rhs.rB_->secondStage() && !rhs.rB_->fixed() ) {
200  rB_->secondStage(kTRUE);
201  rB_->initValue(0.0);
202  }
203  if ( rhs.deltaB_->secondStage() && !rhs.deltaB_->fixed() ) {
204  deltaB_->secondStage(kTRUE);
205  deltaB_->initValue(0.0);
206  }
207  }
209  if ( useGlobalADSPars_ ) {
212  } else {
213  rD_ = new LauParameter("rD", rhs.rD_->value(), minMagnitude_, maxMagnitude_, rhs.rD_->fixed());
214  if ( rhs.rD_->blind() ) {
215  const LauBlind* blinder = rhs.rD_->blinder();
216  rD_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
217  }
218  deltaD_ = new LauParameter("deltaD", rhs.deltaD_->value(), minPhase_, maxPhase_, rhs.deltaD_->fixed());
219  if ( rhs.deltaD_->blind() ) {
220  const LauBlind* blinder = rhs.deltaD_->blinder();
221  deltaD_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
222  }
223  if ( rhs.rD_->secondStage() && !rhs.rD_->fixed() ) {
224  rD_->secondStage(kTRUE);
225  rD_->initValue(0.0);
226  }
227  if ( rhs.deltaD_->secondStage() && !rhs.deltaD_->fixed() ) {
228  deltaD_->secondStage(kTRUE);
229  deltaD_->initValue(0.0);
230  }
231  }
232  }
233  }
234 }
235 
236 void LauPolarGammaCPCoeffSet::adjustName(LauParameter* par, const TString& oldBaseName)
237 {
238  if ( ( par == gamma_ && useGlobalGamma_ ) ||
239  ( par == rD_ && useGlobalADSPars_ ) ||
240  ( par == deltaD_ && useGlobalADSPars_ ) ) {
241  // for global parameters we do not want to adjust their names
242  return;
243  } else {
244  LauAbsCoeffSet::adjustName(par,oldBaseName);
245  }
246 }
247 
248 std::vector<LauParameter*> LauPolarGammaCPCoeffSet::getParameters()
249 {
250  std::vector<LauParameter*> pars;
251  pars.push_back(x_);
252  pars.push_back(y_);
253  if ( !gamma_->fixed() ) {
254  pars.push_back(gamma_);
255  }
257  if ( !rB_->fixed() ) {
258  pars.push_back(rB_);
259  }
260  if ( !deltaB_->fixed() ) {
261  pars.push_back(deltaB_);
262  }
263  }
265  if ( !rD_->fixed() ) {
266  pars.push_back(rD_);
267  }
268  if ( !deltaD_->fixed() ) {
269  pars.push_back(deltaD_);
270  }
271  }
272  return pars;
273 }
274 
276 {
277  std::cout<<"INFO in LauPolarGammaCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
278  std::cout<<"x = "<<x_->value()<<",\t";
279  std::cout<<"y = "<<y_->value()<<",\t";
281  std::cout<<"rB = "<<rB_->value()<<",\t";
282  std::cout<<"deltaB = "<<deltaB_->value()<<",\t";
283  }
285  std::cout<<"rD = "<<rD_->value()<<",\t";
286  std::cout<<"deltaD = "<<deltaD_->value()<<",\t";
287  }
288  std::cout<<"gamma = "<<gamma_->value()<<"."<<std::endl;
289 }
290 
291 void LauPolarGammaCPCoeffSet::printTableHeading(std::ostream& stream) const
292 {
293  switch ( decayType_ ) {
294  case GLW_CPOdd :
295  stream<<"\\begin{tabular}{|l|c|c|c|c|c|}"<<std::endl;
296  stream<<"\\hline"<<std::endl;
297  stream<<"Component & Real Part & Imaginary Part & $r_B$ & $\\delta_B$ & $\\gamma$ \\\\"<<std::endl;
298  break;
299  case GLW_CPEven :
300  stream<<"\\begin{tabular}{|l|c|c|c|c|c|}"<<std::endl;
301  stream<<"\\hline"<<std::endl;
302  stream<<"Component & Real Part & Imaginary Part & $r_B$ & $\\delta_B$ & $\\gamma$ \\\\"<<std::endl;
303  break;
304  case ADS_Favoured :
305  stream<<"\\begin{tabular}{|l|c|c|c|c|c|c|c|}"<<std::endl;
306  stream<<"\\hline"<<std::endl;
307  stream<<"Component & Real Part & Imaginary Part & $r_B$ & $\\delta_B$ & $r_D$ & $\\delta_D$ & $\\gamma$ \\\\"<<std::endl;
308  break;
309  case ADS_Suppressed :
310  stream<<"\\begin{tabular}{|l|c|c|c|c|c|c|c|}"<<std::endl;
311  stream<<"\\hline"<<std::endl;
312  stream<<"Component & Real Part & Imaginary Part & $r_B$ & $\\delta_B$ & $r_D$ & $\\delta_D$ & $\\gamma$ \\\\"<<std::endl;
313  break;
314  case GLW_CPOdd_btouOnly :
315  stream<<"\\begin{tabular}{|l|c|c|c|}"<<std::endl;
316  stream<<"\\hline"<<std::endl;
317  stream<<"Component & Real Part & Imaginary Part & $\\gamma$ \\\\"<<std::endl;
318  break;
319  case GLW_CPEven_btouOnly :
320  stream<<"\\begin{tabular}{|l|c|c|c|}"<<std::endl;
321  stream<<"\\hline"<<std::endl;
322  stream<<"Component & Real Part & Imaginary Part & $\\gamma$ \\\\"<<std::endl;
323  break;
324  case ADS_Favoured_btouOnly :
325  stream<<"\\begin{tabular}{|l|c|c|c|c|c|}"<<std::endl;
326  stream<<"\\hline"<<std::endl;
327  stream<<"Component & Real Part & Imaginary Part & $r_D$ & $\\delta_D$ & $\\gamma$ \\\\"<<std::endl;
328  break;
330  stream<<"\\begin{tabular}{|l|c|c|c|}"<<std::endl;
331  stream<<"\\hline"<<std::endl;
332  stream<<"Component & Real Part & Imaginary Part & $\\gamma$ \\\\"<<std::endl;
333  break;
334  }
335  stream<<"\\hline"<<std::endl;
336 }
337 
338 void LauPolarGammaCPCoeffSet::printTableRow(std::ostream& stream) const
339 {
340  LauPrint print;
341  TString resName = this->name();
342  resName = resName.ReplaceAll("_", "\\_");
343  stream<<resName<<" & $";
344  print.printFormat(stream, x_->value());
345  stream<<" \\pm ";
346  print.printFormat(stream, x_->error());
347  stream<<"$ & $";
348  print.printFormat(stream, y_->value());
349  stream<<" \\pm ";
350  print.printFormat(stream, y_->error());
351  stream<<"$ & $";
353  print.printFormat(stream, rB_->value());
354  stream<<" \\pm ";
355  print.printFormat(stream, rB_->error());
356  stream<<"$ & $";
357  print.printFormat(stream, deltaB_->value());
358  stream<<" \\pm ";
359  print.printFormat(stream, deltaB_->error());
360  stream<<"$ & $";
361  }
363  print.printFormat(stream, rD_->value());
364  stream<<" \\pm ";
365  print.printFormat(stream, rD_->error());
366  stream<<"$ & $";
367  print.printFormat(stream, deltaD_->value());
368  stream<<" \\pm ";
369  print.printFormat(stream, deltaD_->error());
370  stream<<"$ & $";
371  }
372  print.printFormat(stream, gamma_->value());
373  stream<<" \\pm ";
374  print.printFormat(stream, gamma_->error());
375  stream<<"$ \\\\"<<std::endl;
376 }
377 
379 {
380  if (x_->fixed() == kFALSE) {
381  // Choose a value for "X" between -3.0 and 3.0
382  Double_t value = LauRandom::zeroSeedRandom()->Rndm()*6.0 - 3.0;
383  x_->initValue(value); x_->value(value);
384  }
385  if (y_->fixed() == kFALSE) {
386  // Choose a value for "Y" between -3.0 and 3.0
387  Double_t value = LauRandom::zeroSeedRandom()->Rndm()*6.0 - 3.0;
388  y_->initValue(value); y_->value(value);
389  }
390  if (gamma_->fixed() == kFALSE && gamma_->secondStage() == kFALSE) {
391  // Choose a value for "gamma" between +-pi
393  gamma_->initValue(value); gamma_->value(value);
394  }
396  if (rB_->fixed() == kFALSE && rB_->secondStage() == kFALSE) {
397  // Choose a value for "rB" between 0.0 and 2.0
398  Double_t value = LauRandom::zeroSeedRandom()->Rndm()*2.0;
399  rB_->initValue(value); rB_->value(value);
400  }
401  if (deltaB_->fixed() == kFALSE && deltaB_->secondStage() == kFALSE) {
402  // Choose a value for "deltaB" between +- pi
404  deltaB_->initValue(value); deltaB_->value(value);
405  }
406  }
408  if (rD_->fixed() == kFALSE && rD_->secondStage() == kFALSE) {
409  // Choose a value for "rD" between 0.0 and 2.0
410  Double_t value = LauRandom::zeroSeedRandom()->Rndm()*2.0;
411  rD_->initValue(value); rD_->value(value);
412  }
413  if (deltaD_->fixed() == kFALSE && deltaD_->secondStage() == kFALSE) {
414  // Choose a value for "deltaD" between +- pi
416  deltaD_->initValue(value); deltaD_->value(value);
417  }
418  }
419 }
420 
422 {
423  // retrieve the current values from the parameters
424  Double_t gammaVal = gamma_->value();
425  Double_t rBVal = 0.0;
426  Double_t deltaBVal = 0.0;
427  Double_t genDeltaB = 0.0;
428  Double_t rDVal = 0.0;
429  Double_t deltaDVal = 0.0;
430  Double_t genDeltaD = 0.0;
432  rBVal = rB_->value();
433  deltaBVal = deltaB_->value();
434  genDeltaB = deltaB_->genValue();
435  }
437  rDVal = rD_->value();
438  deltaDVal = deltaD_->value();
439  genDeltaD = deltaD_->genValue();
440  }
441 
442 
443  // Check whether we have a negative magnitude.
444  // If so make it positive and add pi to the phases.
445  if (rBVal < 0.0) {
446  rBVal *= -1.0;
447  deltaBVal += LauConstants::pi;
448  }
449  if (rDVal < 0.0) {
450  rDVal *= -1.0;
451  deltaDVal += LauConstants::pi;
452  }
453 
454  // Check now whether the phases lie in the right range (-pi to pi).
455  Bool_t deltaBWithinRange(kFALSE);
456  Bool_t deltaDWithinRange(kFALSE);
457  Bool_t gammaWithinRange(kFALSE);
458  while ( deltaBWithinRange == kFALSE ) {
459  if (deltaBVal > -LauConstants::pi && deltaBVal <= LauConstants::pi) {
460  deltaBWithinRange = kTRUE;
461  } else {
462  // Not within the specified range
463  if (deltaBVal > LauConstants::pi) {
464  deltaBVal -= LauConstants::twoPi;
465  } else if (deltaBVal <= -LauConstants::pi) {
466  deltaBVal += LauConstants::twoPi;
467  }
468  }
469  }
470 
471  while ( deltaDWithinRange == kFALSE ) {
472  if (deltaDVal > -LauConstants::pi && deltaDVal <= LauConstants::pi) {
473  deltaDWithinRange = kTRUE;
474  } else {
475  // Not within the specified range
476  if (deltaDVal > LauConstants::pi) {
477  deltaDVal -= LauConstants::twoPi;
478  } else if (deltaDVal <= -LauConstants::pi) {
479  deltaDVal += LauConstants::twoPi;
480  }
481  }
482  }
483 
484  while ( gammaWithinRange == kFALSE ) {
485  if (gammaVal > -LauConstants::pi && gammaVal <= LauConstants::pi) {
486  gammaWithinRange = kTRUE;
487  } else {
488  // Not within the specified range
489  if (gammaVal > LauConstants::pi) {
490  gammaVal -= LauConstants::twoPi;
491  } else if (gammaVal <= -LauConstants::pi) {
492  gammaVal += LauConstants::twoPi;
493  }
494  }
495  }
496 
497  // To resolve the two-fold ambiguity in gamma and deltaB we require gamma to be in the range 0-pi
499  if (gammaVal < 0.0) {
500  if (deltaBVal <= 0.0) {
501  gammaVal += LauConstants::pi;
502  deltaBVal += LauConstants::pi;
503  } else {
504  gammaVal += LauConstants::pi;
505  deltaBVal -= LauConstants::pi;
506  }
507  }
508  }
509 
510  // A further problem can occur when the generated phase is close to -pi or pi.
511  // The phase can wrap over to the other end of the scale -
512  // this leads to artificially large pulls so we wrap it back.
513  Double_t diff = deltaBVal - genDeltaB;
514  if (diff > LauConstants::pi) {
515  deltaBVal -= LauConstants::twoPi;
516  } else if (diff < -LauConstants::pi) {
517  deltaBVal += LauConstants::twoPi;
518  }
519 
520  diff = deltaDVal - genDeltaD;
521  if (diff > LauConstants::pi) {
522  deltaDVal -= LauConstants::twoPi;
523  } else if (diff < -LauConstants::pi) {
524  deltaDVal += LauConstants::twoPi;
525  }
526 
527  // finally store the new values in the parameters
528  // and update the pulls
529  gamma_->value(gammaVal);
530  gamma_->updatePull();
532  rB_->value(rBVal);
533  rB_->updatePull();
534  deltaB_->value(deltaBVal);
535  deltaB_->updatePull();
536  }
538  rD_->value(rDVal);
539  rD_->updatePull();
540  deltaD_->value(deltaDVal);
541  deltaD_->updatePull();
542  }
543 }
544 
546 {
547  this->updateAmplitudes();
548  return particleCoeff_;
549 }
550 
552 {
553  this->updateAmplitudes();
554  return antiparticleCoeff_;
555 }
556 
558 {
560 
561  const Double_t gammaVal = gamma_->unblindValue();
562 
563  switch ( decayType_ ) {
564 
565  case GLW_CPOdd :
566  {
567  const Double_t rBVal = rB_->unblindValue();
568  const Double_t deltaBVal = deltaB_->unblindValue();
569  cpPart_.setRealImagPart( 1.0 - rBVal*TMath::Cos(deltaBVal + gammaVal), -rBVal*TMath::Sin(deltaBVal + gammaVal) );
570  cpAntiPart_.setRealImagPart( 1.0 - rBVal*TMath::Cos(deltaBVal - gammaVal), -rBVal*TMath::Sin(deltaBVal - gammaVal) );
571  break;
572  }
573 
574  case GLW_CPEven :
575  {
576  const Double_t rBVal = rB_->unblindValue();
577  const Double_t deltaBVal = deltaB_->unblindValue();
578  cpPart_.setRealImagPart( 1.0 + rBVal*TMath::Cos(deltaBVal + gammaVal), rBVal*TMath::Sin(deltaBVal + gammaVal) );
579  cpAntiPart_.setRealImagPart( 1.0 + rBVal*TMath::Cos(deltaBVal - gammaVal), rBVal*TMath::Sin(deltaBVal - gammaVal) );
580  break;
581  }
582 
583  case ADS_Favoured :
584  {
585  const Double_t rBVal = rB_->unblindValue();
586  const Double_t deltaBVal = deltaB_->unblindValue();
587  const Double_t rDVal = rD_->unblindValue();
588  const Double_t deltaDVal = deltaD_->unblindValue();
589  cpPart_.setRealImagPart( 1.0 + rBVal*rDVal*TMath::Cos(deltaBVal - deltaDVal + gammaVal), rBVal*rDVal*TMath::Sin(deltaBVal - deltaDVal + gammaVal) );
590  cpAntiPart_.setRealImagPart( 1.0 + rBVal*rDVal*TMath::Cos(deltaBVal - deltaDVal - gammaVal), rBVal*rDVal*TMath::Sin(deltaBVal - deltaDVal - gammaVal) );
591  break;
592  }
593 
594  case ADS_Suppressed :
595  {
596  const Double_t rBVal = rB_->unblindValue();
597  const Double_t deltaBVal = deltaB_->unblindValue();
598  const Double_t rDVal = rD_->unblindValue();
599  const Double_t deltaDVal = deltaD_->unblindValue();
600  cpPart_.setRealImagPart( rDVal*TMath::Cos(-deltaDVal) + rBVal*TMath::Cos(deltaBVal + gammaVal), rDVal*TMath::Sin(-deltaDVal) + rBVal*TMath::Sin(deltaBVal + gammaVal) );
601  cpAntiPart_.setRealImagPart( rDVal*TMath::Cos(-deltaDVal) + rBVal*TMath::Cos(deltaBVal - gammaVal), rDVal*TMath::Sin(-deltaDVal) + rBVal*TMath::Sin(deltaBVal - gammaVal) );
602  break;
603  }
604 
605  case GLW_CPOdd_btouOnly :
606  nonCPPart_.rescale(-1.0);
607  cpPart_.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
608  cpAntiPart_.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
609  break;
610 
611  case GLW_CPEven_btouOnly :
612  cpPart_.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
613  cpAntiPart_.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
614  break;
615 
616  case ADS_Favoured_btouOnly :
617  {
618  const Double_t rDVal = rD_->unblindValue();
619  const Double_t deltaDVal = deltaD_->unblindValue();
620  cpPart_.setRealImagPart( rDVal * TMath::Cos( -deltaDVal + gammaVal ), rDVal * TMath::Sin( -deltaDVal + gammaVal ) );
621  cpAntiPart_.setRealImagPart( rDVal * TMath::Cos( -deltaDVal - gammaVal ), rDVal * TMath::Sin( -deltaDVal - gammaVal ) );
622  break;
623  }
624 
626  cpPart_.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
627  cpAntiPart_.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
628  break;
629 
630  }
631 
634 }
635 
637 {
638  std::cerr << "ERROR in LauPolarGammaCPCoeffSet::setCoeffValues : Method not supported by this class - too many parameters" << std::endl;
639 }
640 
642 {
643  // set the name
644  TString parName(this->baseName()); parName += "_ACP";
645  acp_.name(parName);
646 
647  // work out the ACP value
648  LauComplex nonCPPart( x_->value(), y_->value() );
649  LauComplex cpPart;
650  LauComplex cpAntiPart;
651 
652  const Double_t gammaVal = gamma_->value();
653 
654  switch ( decayType_ ) {
655 
656  case GLW_CPOdd :
657  {
658  const Double_t rBVal = rB_->value();
659  const Double_t deltaBVal = deltaB_->value();
660  cpPart.setRealImagPart( 1.0 - rBVal*TMath::Cos(deltaBVal + gammaVal), -rBVal*TMath::Sin(deltaBVal + gammaVal) );
661  cpAntiPart.setRealImagPart( 1.0 - rBVal*TMath::Cos(deltaBVal - gammaVal), -rBVal*TMath::Sin(deltaBVal - gammaVal) );
662  break;
663  }
664 
665  case GLW_CPEven :
666  {
667  const Double_t rBVal = rB_->value();
668  const Double_t deltaBVal = deltaB_->value();
669  cpPart.setRealImagPart( 1.0 + rBVal*TMath::Cos(deltaBVal + gammaVal), rBVal*TMath::Sin(deltaBVal + gammaVal) );
670  cpAntiPart.setRealImagPart( 1.0 + rBVal*TMath::Cos(deltaBVal - gammaVal), rBVal*TMath::Sin(deltaBVal - gammaVal) );
671  break;
672  }
673 
674  case ADS_Favoured :
675  {
676  const Double_t rBVal = rB_->value();
677  const Double_t deltaBVal = deltaB_->value();
678  const Double_t rDVal = rD_->value();
679  const Double_t deltaDVal = deltaD_->value();
680  cpPart.setRealImagPart( 1.0 + rBVal*rDVal*TMath::Cos(deltaBVal - deltaDVal + gammaVal), rBVal*rDVal*TMath::Sin(deltaBVal - deltaDVal + gammaVal) );
681  cpAntiPart.setRealImagPart( 1.0 + rBVal*rDVal*TMath::Cos(deltaBVal - deltaDVal - gammaVal), rBVal*rDVal*TMath::Sin(deltaBVal - deltaDVal - gammaVal) );
682  break;
683  }
684 
685  case ADS_Suppressed :
686  {
687  const Double_t rBVal = rB_->value();
688  const Double_t deltaBVal = deltaB_->value();
689  const Double_t rDVal = rD_->value();
690  const Double_t deltaDVal = deltaD_->value();
691  cpPart.setRealImagPart( rDVal*TMath::Cos(-deltaDVal) + rBVal*TMath::Cos(deltaBVal + gammaVal), rDVal*TMath::Sin(-deltaDVal) + rBVal*TMath::Sin(deltaBVal + gammaVal) );
692  cpAntiPart.setRealImagPart( rDVal*TMath::Cos(-deltaDVal) + rBVal*TMath::Cos(deltaBVal - gammaVal), rDVal*TMath::Sin(-deltaDVal) + rBVal*TMath::Sin(deltaBVal - gammaVal) );
693  break;
694  }
695 
696  case GLW_CPOdd_btouOnly :
697  nonCPPart.rescale(-1.0);
698  cpPart.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
699  cpAntiPart.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
700  break;
701 
702  case GLW_CPEven_btouOnly :
703  cpPart.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
704  cpAntiPart.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
705  break;
706 
707  case ADS_Favoured_btouOnly :
708  {
709  const Double_t rDVal = rD_->value();
710  const Double_t deltaDVal = deltaD_->value();
711  cpPart.setRealImagPart( rDVal * TMath::Cos( -deltaDVal + gammaVal ), rDVal * TMath::Sin( -deltaDVal + gammaVal ) );
712  cpAntiPart.setRealImagPart( rDVal * TMath::Cos( -deltaDVal - gammaVal ), rDVal * TMath::Sin( -deltaDVal - gammaVal ) );
713  break;
714  }
715 
717  cpPart.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
718  cpAntiPart.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
719  break;
720 
721 
722  }
723 
724  const LauComplex partCoeff = nonCPPart * cpPart;
725  const LauComplex antiCoeff = nonCPPart * cpAntiPart;
726 
727  const Double_t numer = antiCoeff.abs2() - partCoeff.abs2();
728  const Double_t denom = antiCoeff.abs2() + partCoeff.abs2();
729  const Double_t value = numer/denom;
730 
731  // is it fixed?
732  const Bool_t fixed = gamma_->fixed();
733  acp_.fixed(fixed);
734 
735  // we can't work out the error without the covariance matrix
736  const Double_t error(0.0);
737 
738  // set the value and error
739  acp_.valueAndErrors(value,error);
740 
741  return acp_;
742 }
743 
744 LauAbsCoeffSet* LauPolarGammaCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
745 {
746  LauAbsCoeffSet* clone(0);
747  if ( cloneOption == All || cloneOption == TieRealPart || cloneOption == TieImagPart || cloneOption == TieCPPars ) {
748  clone = new LauPolarGammaCPCoeffSet( *this, cloneOption, constFactor );
749  clone->name( newName );
750  } else {
751  std::cerr << "ERROR in LauPolarGammaCPCoeffSet::createClone : Invalid clone option" << std::endl;
752  }
753  return clone;
754 }
755 
static LauParameter * rDGlobal_
the magnitude of the ratio of the favoured and suppressed D-decay amplitudes (shared by multiple reso...
LauParameter * rD_
the magnitude of the ratio of the favoured and suppressed D-decay amplitudes
static Double_t maxPhase_
Maximum allowed value of phase parameters.
Class for defining a complex coefficient useful for extracting the CKM angle gamma from B -&gt; D h h Da...
Bool_t fixed() const
Check whether the parameter is fixed or floated.
LauComplex nonCPPart_
The b -&gt; c part of the complex coefficient.
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)
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
const TString & name() const
The parameter name.
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.
static Double_t maxRealImagPart_
Maximum allowed value of real/imaginary part parameters.
const DecayType decayType_
The type of the D decay.
LauParameter * deltaD_
the relative strong phase of the favoured and suppressed D-decay amplitudes
LauParameter acp_
The CP asymmetry.
const Bool_t useGlobalADSPars_
Whether the global rD and deltaD are used for this resonance.
virtual void adjustName(LauParameter *par, const TString &oldBaseName)
Prepend the base name and index to the name of a parameter.
LauParameter * gamma_
the relative CP-violating (weak) phase of the b -&gt; u and b -&gt; c amplitudes
File containing declaration of LauPrint class.
LauPolarGammaCPCoeffSet(const TString &compName, const DecayType decayType, const Double_t x, const Double_t y, const Double_t rB, const Double_t deltaB, const Double_t gamma, const Double_t rD, const Double_t deltaD, const Bool_t xFixed, const Bool_t yFixed, const Bool_t rBFixed, const Bool_t deltaBFixed, const Bool_t gammaFixed, const Bool_t rDFixed, const Bool_t deltaDFixed, const Bool_t rBSecondStage=kFALSE, const Bool_t deltaBSecondStage=kFALSE, const Bool_t gammaSecondStage=kFALSE, const Bool_t rDSecondStage=kFALSE, const Bool_t deltaDSecondStage=kFALSE, const Bool_t useGlobalGamma=kFALSE, const Bool_t useGlobalADSPars=kFALSE)
Constructor.
LauComplex cpAntiPart_
The b -&gt; u part of the complex coefficient for the antiparticle.
Class to define various output print commands.
Definition: LauPrint.hh:29
const Bool_t useGlobalGamma_
Whether the global gamma is used for this resonance.
virtual void printTableHeading(std::ostream &stream) const
Print the column headings for a results table.
CloneOption
Options for cloning operation.
Double_t abs2() const
Obtain the square of the absolute value of the complex number.
Definition: LauComplex.hh:232
Bool_t clone() const
Check whether is a clone or not.
LauComplex cpPart_
The b -&gt; u part of the complex coefficient for the particle.
Bool_t blind() const
The blinding state.
virtual const LauComplex & antiparticleCoeff()
Retrieve the complex coefficient for an antiparticle.
LauParameter * x_
The real part of the b -&gt; c amplitude.
virtual void randomiseInitValues()
Randomise the starting values of the parameters for a fit.
static LauParameter * deltaDGlobal_
the relative strong phase of the favoured and suppressed D-decay amplitudes (shared by multiple reson...
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...
LauParameter * y_
The imaginary part of the b -&gt; c amplitude.
LauParameter * rB_
the magnitude of the ratio of the b -&gt; u and b -&gt; c amplitudes
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 void printParValues() const
Print the current values of the parameters.
static Double_t minRealImagPart_
Minimum allowed value of real/imaginary part parameters.
File containing LauRandom namespace.
const LauBlind * blinder() const
Access the blinder object.
LauParameter * deltaB_
the relative CP-conserving (strong) phase of the b -&gt; u and b -&gt; c amplitudes
virtual void finaliseValues()
Make sure values are in &quot;standard&quot; ranges, e.g. phases should be between -pi and pi.
void updateAmplitudes()
Update the amplitudes based on the new values of the parameters.
virtual void printTableRow(std::ostream &stream) const
Print the parameters of the complex coefficient as a row in the results table.
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
static LauParameter * gammaGlobal_
The CP-violating phase (shared by multiple resonances)
Double_t initValue() const
The initial value of the parameter.
void blindParameter(const TString &blindingString, const Double_t width)
Blind the parameter.
File containing declaration of LauPolarGammaCPCoeffSet class.
File containing LauConstants namespace.
virtual std::vector< LauParameter * > getParameters()
Retrieve the parameters of the coefficient, e.g. so that they can be loaded into a fit...
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.
virtual LauParameter acp()
Calculate the CP asymmetry.
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.
static Double_t minPhase_
Minimum allowed value of phase parameters.
Double_t value() const
The value of the parameter.
virtual LauAbsCoeffSet * createClone(const TString &newName, CloneOption cloneOption=All, Double_t constFactor=1.0)
Create a clone of the coefficient set.
static Double_t minMagnitude_
Minimum allowed value of magnitude parameters.
Double_t blindingWidth() const
Obtain the Gaussian width.
Definition: LauBlind.hh:68
LauComplex particleCoeff_
The particle complex coefficient.
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.
virtual void adjustName(LauParameter *par, const TString &oldBaseName)
Prepend the base name and index to the name of a parameter.
DecayType
The possible D decay modes.