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