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