laura is hosted by Hepforge, IPPP Durham
Laura++  3.6.0
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 
30 
31 #include "LauComplex.hh"
32 #include "LauConstants.hh"
33 #include "LauParameter.hh"
34 #include "LauPrint.hh"
35 
36 #include "TMath.h"
37 #include "TRandom.h"
38 
39 #include <fstream>
40 #include <iostream>
41 #include <vector>
42 
46 
48  const DecayType decayType,
49  const Double_t x,
50  const Double_t y,
51  const Double_t rB,
52  const Double_t deltaB,
53  const Double_t gamma,
54  const Double_t rD,
55  const Double_t deltaD,
56  const Bool_t xFixed,
57  const Bool_t yFixed,
58  const Bool_t rBFixed,
59  const Bool_t deltaBFixed,
60  const Bool_t gammaFixed,
61  const Bool_t rDFixed,
62  const Bool_t deltaDFixed,
63  const Bool_t rBSecondStage,
64  const Bool_t deltaBSecondStage,
65  const Bool_t gammaSecondStage,
66  const Bool_t rDSecondStage,
67  const Bool_t deltaDSecondStage,
68  const Bool_t useGlobalGamma,
69  const Bool_t useGlobalADSPars ) :
70  LauAbsCoeffSet( compName ),
71  decayType_( decayType ),
72  x_( 0 ),
73  y_( 0 ),
74  rB_( 0 ),
75  deltaB_( 0 ),
76  gamma_( 0 ),
77  rD_( 0 ),
78  deltaD_( 0 ),
79  useGlobalGamma_( useGlobalGamma ),
80  useGlobalADSPars_( useGlobalADSPars ),
81  nonCPPart_( x, y ),
82  cpPart_( 0.0, 0.0 ),
83  cpAntiPart_( 0.0, 0.0 ),
84  particleCoeff_( 0.0, 0.0 ),
85  antiparticleCoeff_( 0.0, 0.0 ),
86  acp_( "ACP", 0.0, -1.0, 1.0 )
87 {
88  // All of the possible D decay types need these two parameters
89  x_ = new LauParameter( "X", x, minRealImagPart_, maxRealImagPart_, xFixed );
90  y_ = new LauParameter( "Y", y, minRealImagPart_, maxRealImagPart_, yFixed );
91 
92  // if we're using a global gamma, create it if it doesn't already exist then set gamma_ to point to it
93  // otherwise create our individual copy of gamma
94  if ( useGlobalGamma_ ) {
95  if ( ! gammaGlobal_ ) {
96  gammaGlobal_ = new LauParameter( "gamma", gamma, minPhase_, maxPhase_, gammaFixed );
98  } else {
100  }
101  } else {
102  gamma_ = new LauParameter( "gamma", gamma, minPhase_, maxPhase_, gammaFixed );
103  }
104  if ( gammaSecondStage && ! gammaFixed ) {
105  gamma_->secondStage( kTRUE );
106  gamma_->initValue( 0.0 );
107  }
108 
109  // which of the other parameter we need depends on the D-decay type
111  decayType_ == GLW_CPEven ) {
112  rB_ = new LauParameter( "rB", rB, minMagnitude_, maxMagnitude_, rBFixed );
113  deltaB_ = new LauParameter( "deltaB", deltaB, minPhase_, maxPhase_, deltaBFixed );
114  if ( rBSecondStage && ! rBFixed ) {
115  rB_->secondStage( kTRUE );
116  rB_->initValue( 0.0 );
117  }
118  if ( deltaBSecondStage && ! deltaBFixed ) {
119  deltaB_->secondStage( kTRUE );
120  deltaB_->initValue( 0.0 );
121  }
122  }
123 
126  if ( useGlobalADSPars_ ) {
127  if ( ! rDGlobal_ ) {
128  rDGlobal_ = new LauParameter( "rD", rD, minMagnitude_, maxMagnitude_, rDFixed );
129  deltaDGlobal_ = new LauParameter( "deltaD", deltaD, minPhase_, maxPhase_, deltaDFixed );
130  rD_ = rDGlobal_;
132  } else {
135  }
136  } else {
137  rD_ = new LauParameter( "rD", rD, minMagnitude_, maxMagnitude_, rDFixed );
138  deltaD_ = new LauParameter( "deltaD", deltaD, minPhase_, maxPhase_, deltaDFixed );
139  }
140  if ( rDSecondStage && ! rDFixed ) {
141  rD_->secondStage( kTRUE );
142  rD_->initValue( 0.0 );
143  }
144  if ( deltaDSecondStage && ! deltaDFixed ) {
145  deltaD_->secondStage( kTRUE );
146  deltaD_->initValue( 0.0 );
147  }
148  }
149 }
150 
152  CloneOption cloneOption,
153  Double_t constFactor ) :
154  LauAbsCoeffSet( rhs.name() ),
155  decayType_( rhs.decayType_ ),
156  x_( 0 ),
157  y_( 0 ),
158  rB_( 0 ),
159  deltaB_( 0 ),
160  gamma_( 0 ),
161  rD_( 0 ),
162  deltaD_( 0 ),
163  useGlobalGamma_( rhs.useGlobalGamma_ ),
164  useGlobalADSPars_( rhs.useGlobalADSPars_ ),
165  nonCPPart_( rhs.nonCPPart_ ),
166  cpPart_( rhs.cpPart_ ),
167  cpAntiPart_( rhs.cpAntiPart_ ),
168  particleCoeff_( rhs.particleCoeff_ ),
169  antiparticleCoeff_( rhs.antiparticleCoeff_ ),
170  acp_( rhs.acp_ )
171 {
172  if ( cloneOption == All || cloneOption == TieRealPart ) {
173  x_ = rhs.x_->createClone( constFactor );
174  } else {
175  x_ = new LauParameter( "X",
176  rhs.x_->value(),
179  rhs.x_->fixed() );
180  if ( rhs.x_->blind() ) {
181  const LauBlind* blinder = rhs.x_->blinder();
183  }
184  }
185 
186  if ( cloneOption == All || cloneOption == TieImagPart ) {
187  y_ = rhs.y_->createClone( constFactor );
188  } else {
189  y_ = new LauParameter( "Y",
190  rhs.y_->value(),
193  rhs.y_->fixed() );
194  if ( rhs.y_->blind() ) {
195  const LauBlind* blinder = rhs.y_->blinder();
197  }
198  }
199 
200  if ( cloneOption == All || cloneOption == TieCPPars ) {
201  gamma_ = rhs.gamma_->createClone( constFactor );
204  rB_ = rhs.rB_->createClone( constFactor );
205  deltaB_ = rhs.deltaB_->createClone( constFactor );
206  }
209  rD_ = rhs.rD_->createClone( constFactor );
210  deltaD_ = rhs.deltaD_->createClone( constFactor );
211  }
212  } else {
213  if ( useGlobalGamma_ ) {
215  } else {
216  gamma_ = new LauParameter( "gamma",
217  rhs.gamma_->value(),
218  minPhase_,
219  maxPhase_,
220  rhs.gamma_->fixed() );
221  if ( rhs.gamma_->blind() ) {
222  const LauBlind* blinder = rhs.gamma_->blinder();
224  }
225  if ( rhs.gamma_->secondStage() && ! rhs.gamma_->fixed() ) {
226  gamma_->secondStage( kTRUE );
227  gamma_->initValue( 0.0 );
228  }
229  }
232  rB_ = new LauParameter( "rB",
233  rhs.rB_->value(),
236  rhs.rB_->fixed() );
237  if ( rhs.rB_->blind() ) {
238  const LauBlind* blinder = rhs.rB_->blinder();
240  }
241  deltaB_ = new LauParameter( "deltaB",
242  rhs.deltaB_->value(),
243  minPhase_,
244  maxPhase_,
245  rhs.deltaB_->fixed() );
246  if ( rhs.deltaB_->blind() ) {
247  const LauBlind* blinder = rhs.deltaB_->blinder();
249  }
250  if ( rhs.rB_->secondStage() && ! rhs.rB_->fixed() ) {
251  rB_->secondStage( kTRUE );
252  rB_->initValue( 0.0 );
253  }
254  if ( rhs.deltaB_->secondStage() && ! rhs.deltaB_->fixed() ) {
255  deltaB_->secondStage( kTRUE );
256  deltaB_->initValue( 0.0 );
257  }
258  }
261  if ( useGlobalADSPars_ ) {
264  } else {
265  rD_ = new LauParameter( "rD",
266  rhs.rD_->value(),
269  rhs.rD_->fixed() );
270  if ( rhs.rD_->blind() ) {
271  const LauBlind* blinder = rhs.rD_->blinder();
273  }
274  deltaD_ = new LauParameter( "deltaD",
275  rhs.deltaD_->value(),
276  minPhase_,
277  maxPhase_,
278  rhs.deltaD_->fixed() );
279  if ( rhs.deltaD_->blind() ) {
280  const LauBlind* blinder = rhs.deltaD_->blinder();
282  }
283  if ( rhs.rD_->secondStage() && ! rhs.rD_->fixed() ) {
284  rD_->secondStage( kTRUE );
285  rD_->initValue( 0.0 );
286  }
287  if ( rhs.deltaD_->secondStage() && ! rhs.deltaD_->fixed() ) {
288  deltaD_->secondStage( kTRUE );
289  deltaD_->initValue( 0.0 );
290  }
291  }
292  }
293  }
294 }
295 
296 void LauPolarGammaCPCoeffSet::adjustName( LauParameter* par, const TString& oldBaseName )
297 {
298  if ( ( par == gamma_ && useGlobalGamma_ ) || ( par == rD_ && useGlobalADSPars_ ) ||
299  ( par == deltaD_ && useGlobalADSPars_ ) ) {
300  // for global parameters we do not want to adjust their names
301  return;
302  } else {
303  LauAbsCoeffSet::adjustName( par, oldBaseName );
304  }
305 }
306 
307 std::vector<LauParameter*> LauPolarGammaCPCoeffSet::getParameters()
308 {
309  std::vector<LauParameter*> pars;
310  pars.push_back( x_ );
311  pars.push_back( y_ );
312  if ( ! gamma_->fixed() ) {
313  pars.push_back( gamma_ );
314  }
316  decayType_ == GLW_CPEven ) {
317  if ( ! rB_->fixed() ) {
318  pars.push_back( rB_ );
319  }
320  if ( ! deltaB_->fixed() ) {
321  pars.push_back( deltaB_ );
322  }
323  }
326  if ( ! rD_->fixed() ) {
327  pars.push_back( rD_ );
328  }
329  if ( ! deltaD_->fixed() ) {
330  pars.push_back( deltaD_ );
331  }
332  }
333  return pars;
334 }
335 
337 {
338  std::cout << "INFO in LauPolarGammaCPCoeffSet::printParValues : Component \"" << this->name()
339  << "\" has ";
340  std::cout << "x = " << x_->value() << ",\t";
341  std::cout << "y = " << y_->value() << ",\t";
343  decayType_ == GLW_CPEven ) {
344  std::cout << "rB = " << rB_->value() << ",\t";
345  std::cout << "deltaB = " << deltaB_->value() << ",\t";
346  }
349  std::cout << "rD = " << rD_->value() << ",\t";
350  std::cout << "deltaD = " << deltaD_->value() << ",\t";
351  }
352  std::cout << "gamma = " << gamma_->value() << "." << std::endl;
353 }
354 
355 void LauPolarGammaCPCoeffSet::printTableHeading( std::ostream& stream ) const
356 {
357  switch ( decayType_ ) {
358  case GLW_CPOdd :
359  stream << "\\begin{tabular}{|l|c|c|c|c|c|}" << std::endl;
360  stream << "\\hline" << std::endl;
361  stream << "Component & Real Part & Imaginary Part & $r_B$ & $\\delta_B$ & $\\gamma$ \\\\"
362  << std::endl;
363  break;
364  case GLW_CPEven :
365  stream << "\\begin{tabular}{|l|c|c|c|c|c|}" << std::endl;
366  stream << "\\hline" << std::endl;
367  stream << "Component & Real Part & Imaginary Part & $r_B$ & $\\delta_B$ & $\\gamma$ \\\\"
368  << std::endl;
369  break;
370  case ADS_Favoured :
371  stream << "\\begin{tabular}{|l|c|c|c|c|c|c|c|}" << std::endl;
372  stream << "\\hline" << std::endl;
373  stream << "Component & Real Part & Imaginary Part & $r_B$ & $\\delta_B$ & $r_D$ & $\\delta_D$ & $\\gamma$ \\\\"
374  << std::endl;
375  break;
376  case ADS_Suppressed :
377  stream << "\\begin{tabular}{|l|c|c|c|c|c|c|c|}" << std::endl;
378  stream << "\\hline" << std::endl;
379  stream << "Component & Real Part & Imaginary Part & $r_B$ & $\\delta_B$ & $r_D$ & $\\delta_D$ & $\\gamma$ \\\\"
380  << std::endl;
381  break;
382  case GLW_CPOdd_btouOnly :
383  stream << "\\begin{tabular}{|l|c|c|c|}" << std::endl;
384  stream << "\\hline" << std::endl;
385  stream << "Component & Real Part & Imaginary Part & $\\gamma$ \\\\" << std::endl;
386  break;
387  case GLW_CPEven_btouOnly :
388  stream << "\\begin{tabular}{|l|c|c|c|}" << std::endl;
389  stream << "\\hline" << std::endl;
390  stream << "Component & Real Part & Imaginary Part & $\\gamma$ \\\\" << std::endl;
391  break;
392  case ADS_Favoured_btouOnly :
393  stream << "\\begin{tabular}{|l|c|c|c|c|c|}" << std::endl;
394  stream << "\\hline" << std::endl;
395  stream << "Component & Real Part & Imaginary Part & $r_D$ & $\\delta_D$ & $\\gamma$ \\\\"
396  << std::endl;
397  break;
399  stream << "\\begin{tabular}{|l|c|c|c|}" << std::endl;
400  stream << "\\hline" << std::endl;
401  stream << "Component & Real Part & Imaginary Part & $\\gamma$ \\\\" << std::endl;
402  break;
403  }
404  stream << "\\hline" << std::endl;
405 }
406 
407 void LauPolarGammaCPCoeffSet::printTableRow( std::ostream& stream ) const
408 {
409  LauPrint print;
410  TString resName = this->name();
411  resName = resName.ReplaceAll( "_", "\\_" );
412  stream << resName << " & $";
413  print.printFormat( stream, x_->value() );
414  stream << " \\pm ";
415  print.printFormat( stream, x_->error() );
416  stream << "$ & $";
417  print.printFormat( stream, y_->value() );
418  stream << " \\pm ";
419  print.printFormat( stream, y_->error() );
420  stream << "$ & $";
422  decayType_ == GLW_CPEven ) {
423  print.printFormat( stream, rB_->value() );
424  stream << " \\pm ";
425  print.printFormat( stream, rB_->error() );
426  stream << "$ & $";
427  print.printFormat( stream, deltaB_->value() );
428  stream << " \\pm ";
429  print.printFormat( stream, deltaB_->error() );
430  stream << "$ & $";
431  }
434  print.printFormat( stream, rD_->value() );
435  stream << " \\pm ";
436  print.printFormat( stream, rD_->error() );
437  stream << "$ & $";
438  print.printFormat( stream, deltaD_->value() );
439  stream << " \\pm ";
440  print.printFormat( stream, deltaD_->error() );
441  stream << "$ & $";
442  }
443  print.printFormat( stream, gamma_->value() );
444  stream << " \\pm ";
445  print.printFormat( stream, gamma_->error() );
446  stream << "$ \\\\" << std::endl;
447 }
448 
450 {
451  if ( x_->fixed() == kFALSE ) {
452  // Choose a value for "X" between -3.0 and 3.0
453  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm() * 6.0 - 3.0;
454  x_->initValue( value );
455  x_->value( value );
456  }
457  if ( y_->fixed() == kFALSE ) {
458  // Choose a value for "Y" between -3.0 and 3.0
459  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm() * 6.0 - 3.0;
460  y_->initValue( value );
461  y_->value( value );
462  }
463  if ( gamma_->fixed() == kFALSE && gamma_->secondStage() == kFALSE ) {
464  // Choose a value for "gamma" between +-pi
467  gamma_->initValue( value );
468  gamma_->value( value );
469  }
471  decayType_ == GLW_CPEven ) {
472  if ( rB_->fixed() == kFALSE && rB_->secondStage() == kFALSE ) {
473  // Choose a value for "rB" between 0.0 and 2.0
474  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm() * 2.0;
475  rB_->initValue( value );
476  rB_->value( value );
477  }
478  if ( deltaB_->fixed() == kFALSE && deltaB_->secondStage() == kFALSE ) {
479  // Choose a value for "deltaB" between +- pi
482  deltaB_->initValue( value );
483  deltaB_->value( value );
484  }
485  }
488  if ( rD_->fixed() == kFALSE && rD_->secondStage() == kFALSE ) {
489  // Choose a value for "rD" between 0.0 and 2.0
490  Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm() * 2.0;
491  rD_->initValue( value );
492  rD_->value( value );
493  }
494  if ( deltaD_->fixed() == kFALSE && deltaD_->secondStage() == kFALSE ) {
495  // Choose a value for "deltaD" between +- pi
498  deltaD_->initValue( value );
499  deltaD_->value( value );
500  }
501  }
502 }
503 
505 {
506  // retrieve the current values from the parameters
507  Double_t gammaVal = gamma_->value();
508  Double_t rBVal = 0.0;
509  Double_t deltaBVal = 0.0;
510  Double_t genDeltaB = 0.0;
511  Double_t rDVal = 0.0;
512  Double_t deltaDVal = 0.0;
513  Double_t genDeltaD = 0.0;
515  decayType_ == GLW_CPEven ) {
516  rBVal = rB_->value();
517  deltaBVal = deltaB_->value();
518  genDeltaB = deltaB_->genValue();
519  }
522  rDVal = rD_->value();
523  deltaDVal = deltaD_->value();
524  genDeltaD = deltaD_->genValue();
525  }
526 
527  // Check whether we have a negative magnitude.
528  // If so make it positive and add pi to the phases.
529  if ( rBVal < 0.0 ) {
530  rBVal *= -1.0;
531  deltaBVal += LauConstants::pi;
532  }
533  if ( rDVal < 0.0 ) {
534  rDVal *= -1.0;
535  deltaDVal += LauConstants::pi;
536  }
537 
538  // Check now whether the phases lie in the right range (-pi to pi).
539  Bool_t deltaBWithinRange( kFALSE );
540  Bool_t deltaDWithinRange( kFALSE );
541  Bool_t gammaWithinRange( kFALSE );
542  while ( deltaBWithinRange == kFALSE ) {
543  if ( deltaBVal > -LauConstants::pi && deltaBVal <= LauConstants::pi ) {
544  deltaBWithinRange = kTRUE;
545  } else {
546  // Not within the specified range
547  if ( deltaBVal > LauConstants::pi ) {
548  deltaBVal -= LauConstants::twoPi;
549  } else if ( deltaBVal <= -LauConstants::pi ) {
550  deltaBVal += LauConstants::twoPi;
551  }
552  }
553  }
554 
555  while ( deltaDWithinRange == kFALSE ) {
556  if ( deltaDVal > -LauConstants::pi && deltaDVal <= LauConstants::pi ) {
557  deltaDWithinRange = kTRUE;
558  } else {
559  // Not within the specified range
560  if ( deltaDVal > LauConstants::pi ) {
561  deltaDVal -= LauConstants::twoPi;
562  } else if ( deltaDVal <= -LauConstants::pi ) {
563  deltaDVal += LauConstants::twoPi;
564  }
565  }
566  }
567 
568  while ( gammaWithinRange == kFALSE ) {
569  if ( gammaVal > -LauConstants::pi && gammaVal <= LauConstants::pi ) {
570  gammaWithinRange = kTRUE;
571  } else {
572  // Not within the specified range
573  if ( gammaVal > LauConstants::pi ) {
574  gammaVal -= LauConstants::twoPi;
575  } else if ( gammaVal <= -LauConstants::pi ) {
576  gammaVal += LauConstants::twoPi;
577  }
578  }
579  }
580 
581  // To resolve the two-fold ambiguity in gamma and deltaB we require gamma to be in the range 0-pi
583  decayType_ == GLW_CPEven ) {
584  if ( gammaVal < 0.0 ) {
585  if ( deltaBVal <= 0.0 ) {
586  gammaVal += LauConstants::pi;
587  deltaBVal += LauConstants::pi;
588  } else {
589  gammaVal += LauConstants::pi;
590  deltaBVal -= LauConstants::pi;
591  }
592  }
593  }
594 
595  // A further problem can occur when the generated phase is close to -pi or pi.
596  // The phase can wrap over to the other end of the scale -
597  // this leads to artificially large pulls so we wrap it back.
598  Double_t diff = deltaBVal - genDeltaB;
599  if ( diff > LauConstants::pi ) {
600  deltaBVal -= LauConstants::twoPi;
601  } else if ( diff < -LauConstants::pi ) {
602  deltaBVal += LauConstants::twoPi;
603  }
604 
605  diff = deltaDVal - genDeltaD;
606  if ( diff > LauConstants::pi ) {
607  deltaDVal -= LauConstants::twoPi;
608  } else if ( diff < -LauConstants::pi ) {
609  deltaDVal += LauConstants::twoPi;
610  }
611 
612  // finally store the new values in the parameters
613  // and update the pulls
614  gamma_->value( gammaVal );
615  gamma_->updatePull();
617  decayType_ == GLW_CPEven ) {
618  rB_->value( rBVal );
619  rB_->updatePull();
620  deltaB_->value( deltaBVal );
621  deltaB_->updatePull();
622  }
625  rD_->value( rDVal );
626  rD_->updatePull();
627  deltaD_->value( deltaDVal );
628  deltaD_->updatePull();
629  }
630 }
631 
633 {
634  this->updateAmplitudes();
635  return particleCoeff_;
636 }
637 
639 {
640  this->updateAmplitudes();
641  return antiparticleCoeff_;
642 }
643 
645 {
647 
648  const Double_t gammaVal = gamma_->unblindValue();
649 
650  switch ( decayType_ ) {
651 
652  case GLW_CPOdd : {
653  const Double_t rBVal = rB_->unblindValue();
654  const Double_t deltaBVal = deltaB_->unblindValue();
655  cpPart_.setRealImagPart( 1.0 - rBVal * TMath::Cos( deltaBVal + gammaVal ),
656  -rBVal * TMath::Sin( deltaBVal + gammaVal ) );
657  cpAntiPart_.setRealImagPart( 1.0 - rBVal * TMath::Cos( deltaBVal - gammaVal ),
658  -rBVal * TMath::Sin( deltaBVal - gammaVal ) );
659  break;
660  }
661 
662  case GLW_CPEven : {
663  const Double_t rBVal = rB_->unblindValue();
664  const Double_t deltaBVal = deltaB_->unblindValue();
665  cpPart_.setRealImagPart( 1.0 + rBVal * TMath::Cos( deltaBVal + gammaVal ),
666  rBVal * TMath::Sin( deltaBVal + gammaVal ) );
667  cpAntiPart_.setRealImagPart( 1.0 + rBVal * TMath::Cos( deltaBVal - gammaVal ),
668  rBVal * TMath::Sin( deltaBVal - gammaVal ) );
669  break;
670  }
671 
672  case ADS_Favoured : {
673  const Double_t rBVal = rB_->unblindValue();
674  const Double_t deltaBVal = deltaB_->unblindValue();
675  const Double_t rDVal = rD_->unblindValue();
676  const Double_t deltaDVal = deltaD_->unblindValue();
677  cpPart_.setRealImagPart( 1.0 + rBVal * rDVal *
678  TMath::Cos( deltaBVal - deltaDVal + gammaVal ),
679  rBVal * rDVal * TMath::Sin( deltaBVal - deltaDVal + gammaVal ) );
681  1.0 + rBVal * rDVal * TMath::Cos( deltaBVal - deltaDVal - gammaVal ),
682  rBVal * rDVal * TMath::Sin( deltaBVal - deltaDVal - gammaVal ) );
683  break;
684  }
685 
686  case ADS_Suppressed : {
687  const Double_t rBVal = rB_->unblindValue();
688  const Double_t deltaBVal = deltaB_->unblindValue();
689  const Double_t rDVal = rD_->unblindValue();
690  const Double_t deltaDVal = deltaD_->unblindValue();
692  rDVal * TMath::Cos( -deltaDVal ) + rBVal * TMath::Cos( deltaBVal + gammaVal ),
693  rDVal * TMath::Sin( -deltaDVal ) + rBVal * TMath::Sin( deltaBVal + gammaVal ) );
695  rDVal * TMath::Cos( -deltaDVal ) + rBVal * TMath::Cos( deltaBVal - gammaVal ),
696  rDVal * TMath::Sin( -deltaDVal ) + rBVal * TMath::Sin( deltaBVal - gammaVal ) );
697  break;
698  }
699 
700  case GLW_CPOdd_btouOnly :
701  nonCPPart_.rescale( -1.0 );
702  cpPart_.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
703  cpAntiPart_.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ),
704  1.0 * TMath::Sin( -gammaVal ) );
705  break;
706 
707  case GLW_CPEven_btouOnly :
708  cpPart_.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
709  cpAntiPart_.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ),
710  1.0 * TMath::Sin( -gammaVal ) );
711  break;
712 
713  case ADS_Favoured_btouOnly : {
714  const Double_t rDVal = rD_->unblindValue();
715  const Double_t deltaDVal = deltaD_->unblindValue();
716  cpPart_.setRealImagPart( rDVal * TMath::Cos( -deltaDVal + gammaVal ),
717  rDVal * TMath::Sin( -deltaDVal + gammaVal ) );
718  cpAntiPart_.setRealImagPart( rDVal * TMath::Cos( -deltaDVal - gammaVal ),
719  rDVal * TMath::Sin( -deltaDVal - gammaVal ) );
720  break;
721  }
722 
724  cpPart_.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
725  cpAntiPart_.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ),
726  1.0 * TMath::Sin( -gammaVal ) );
727  break;
728  }
729 
732 }
733 
735 {
736  std::cerr << "ERROR in LauPolarGammaCPCoeffSet::setCoeffValues : Method not supported by this class - too many parameters"
737  << std::endl;
738 }
739 
741 {
742  // set the name
743  TString parName( this->baseName() );
744  parName += "_ACP";
745  acp_.name( parName );
746 
747  // work out the ACP value
748  LauComplex nonCPPart( x_->value(), y_->value() );
749  LauComplex cpPart;
750  LauComplex cpAntiPart;
751 
752  const Double_t gammaVal = gamma_->value();
753 
754  switch ( decayType_ ) {
755 
756  case GLW_CPOdd : {
757  const Double_t rBVal = rB_->value();
758  const Double_t deltaBVal = deltaB_->value();
759  cpPart.setRealImagPart( 1.0 - rBVal * TMath::Cos( deltaBVal + gammaVal ),
760  -rBVal * TMath::Sin( deltaBVal + gammaVal ) );
761  cpAntiPart.setRealImagPart( 1.0 - rBVal * TMath::Cos( deltaBVal - gammaVal ),
762  -rBVal * TMath::Sin( deltaBVal - gammaVal ) );
763  break;
764  }
765 
766  case GLW_CPEven : {
767  const Double_t rBVal = rB_->value();
768  const Double_t deltaBVal = deltaB_->value();
769  cpPart.setRealImagPart( 1.0 + rBVal * TMath::Cos( deltaBVal + gammaVal ),
770  rBVal * TMath::Sin( deltaBVal + gammaVal ) );
771  cpAntiPart.setRealImagPart( 1.0 + rBVal * TMath::Cos( deltaBVal - gammaVal ),
772  rBVal * TMath::Sin( deltaBVal - gammaVal ) );
773  break;
774  }
775 
776  case ADS_Favoured : {
777  const Double_t rBVal = rB_->value();
778  const Double_t deltaBVal = deltaB_->value();
779  const Double_t rDVal = rD_->value();
780  const Double_t deltaDVal = deltaD_->value();
781  cpPart.setRealImagPart( 1.0 + rBVal * rDVal *
782  TMath::Cos( deltaBVal - deltaDVal + gammaVal ),
783  rBVal * rDVal * TMath::Sin( deltaBVal - deltaDVal + gammaVal ) );
784  cpAntiPart.setRealImagPart(
785  1.0 + rBVal * rDVal * TMath::Cos( deltaBVal - deltaDVal - gammaVal ),
786  rBVal * rDVal * TMath::Sin( deltaBVal - deltaDVal - gammaVal ) );
787  break;
788  }
789 
790  case ADS_Suppressed : {
791  const Double_t rBVal = rB_->value();
792  const Double_t deltaBVal = deltaB_->value();
793  const Double_t rDVal = rD_->value();
794  const Double_t deltaDVal = deltaD_->value();
795  cpPart.setRealImagPart(
796  rDVal * TMath::Cos( -deltaDVal ) + rBVal * TMath::Cos( deltaBVal + gammaVal ),
797  rDVal * TMath::Sin( -deltaDVal ) + rBVal * TMath::Sin( deltaBVal + gammaVal ) );
798  cpAntiPart.setRealImagPart(
799  rDVal * TMath::Cos( -deltaDVal ) + rBVal * TMath::Cos( deltaBVal - gammaVal ),
800  rDVal * TMath::Sin( -deltaDVal ) + rBVal * TMath::Sin( deltaBVal - gammaVal ) );
801  break;
802  }
803 
804  case GLW_CPOdd_btouOnly :
805  nonCPPart.rescale( -1.0 );
806  cpPart.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
807  cpAntiPart.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
808  break;
809 
810  case GLW_CPEven_btouOnly :
811  cpPart.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
812  cpAntiPart.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
813  break;
814 
815  case ADS_Favoured_btouOnly : {
816  const Double_t rDVal = rD_->value();
817  const Double_t deltaDVal = deltaD_->value();
818  cpPart.setRealImagPart( rDVal * TMath::Cos( -deltaDVal + gammaVal ),
819  rDVal * TMath::Sin( -deltaDVal + gammaVal ) );
820  cpAntiPart.setRealImagPart( rDVal * TMath::Cos( -deltaDVal - gammaVal ),
821  rDVal * TMath::Sin( -deltaDVal - gammaVal ) );
822  break;
823  }
824 
826  cpPart.setRealImagPart( 1.0 * TMath::Cos( gammaVal ), 1.0 * TMath::Sin( gammaVal ) );
827  cpAntiPart.setRealImagPart( 1.0 * TMath::Cos( -gammaVal ), 1.0 * TMath::Sin( -gammaVal ) );
828  break;
829  }
830 
831  const LauComplex partCoeff = nonCPPart * cpPart;
832  const LauComplex antiCoeff = nonCPPart * cpAntiPart;
833 
834  const Double_t numer = antiCoeff.abs2() - partCoeff.abs2();
835  const Double_t denom = antiCoeff.abs2() + partCoeff.abs2();
836  const Double_t value = numer / denom;
837 
838  // is it fixed?
839  const Bool_t fixed = gamma_->fixed();
840  acp_.fixed( fixed );
841 
842  // we can't work out the error without the covariance matrix
843  const Double_t error( 0.0 );
844 
845  // set the value and error
847 
848  return acp_;
849 }
850 
852  CloneOption cloneOption,
853  Double_t constFactor )
854 {
855  LauAbsCoeffSet* clone( 0 );
856  if ( cloneOption == All || cloneOption == TieRealPart || cloneOption == TieImagPart ||
857  cloneOption == TieCPPars ) {
858  clone = new LauPolarGammaCPCoeffSet( *this, cloneOption, constFactor );
859  clone->name( newName );
860  } else {
861  std::cerr << "ERROR in LauPolarGammaCPCoeffSet::createClone : Invalid clone option"
862  << std::endl;
863  }
864  return clone;
865 }
static LauParameter * gammaGlobal_
The CP-violating phase (shared by multiple resonances)
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.
LauComplex particleCoeff_
The particle complex coefficient.
LauComplex cpAntiPart_
The b -> u part of the complex coefficient for the antiparticle.
Class for defining the abstract interface for complex coefficient classes.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:49
Double_t unblindValue() const
The unblinded value of the parameter.
virtual LauParameter acp()
Calculate the CP asymmetry.
LauParameter acp_
The CP asymmetry.
Double_t value() const
The value of the parameter.
virtual void adjustName(LauParameter *par, const TString &oldBaseName)
Prepend the base name and index to the name of a parameter.
virtual LauAbsCoeffSet * createClone(const TString &newName, CloneOption cloneOption=All, Double_t constFactor=1.0)
Create a clone of the coefficient set.
DecayType
The possible D decay modes.
static LauParameter * rDGlobal_
the magnitude of the ratio of the favoured and suppressed D-decay amplitudes (shared by multiple reso...
File containing declaration of LauPolarGammaCPCoeffSet class.
const LauBlind * blinder() const
Access the blinder object.
virtual TString name() const
Retrieve the name of the coefficient set.
virtual void adjustName(LauParameter *par, const TString &oldBaseName)
Prepend the base name and index to the name of a parameter.
const Bool_t useGlobalGamma_
Whether the global gamma is used for this resonance.
LauComplex cpPart_
The b -> u part of the complex coefficient for the particle.
File containing declaration of LauParameter class.
const Double_t pi
Pi.
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 maxPhase_
Maximum allowed value of phase parameters.
LauParameter * y_
The imaginary part of the b -> c amplitude.
virtual void printParValues() const
Print the current values of the parameters.
virtual const LauComplex & antiparticleCoeff()
Retrieve the complex coefficient for an antiparticle.
static Double_t maxMagnitude_
Maximum allowed value of magnitude parameters.
virtual const TString & baseName() const
Retrieve the base name of the coefficient set.
Bool_t secondStage() const
Check whether the parameter should be floated only in the second stage of a two stage fit.
const Bool_t useGlobalADSPars_
Whether the global rD and deltaD are used for this resonance.
Class for defining a complex number.
Definition: LauComplex.hh:61
virtual void setCoeffValues(const LauComplex &coeff, const LauComplex &coeffBar, Bool_t init)
Set the parameters based on the complex coefficients for particles and antiparticles.
Bool_t blind() const
The blinding state.
const DecayType decayType_
The type of the D decay.
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.
Class for blinding and unblinding a number based on a blinding string.
Definition: LauBlind.hh:42
void rescale(Double_t scaleVal)
Scale this by a factor.
Definition: LauComplex.hh:301
virtual void printTableRow(std::ostream &stream) const
Print the parameters of the complex coefficient as a row in the results table.
virtual void printTableHeading(std::ostream &stream) const
Print the column headings for a results table.
void setRealImagPart(Double_t realpart, Double_t imagpart)
Set both real and imaginary part.
Definition: LauComplex.hh:324
LauParameter * deltaD_
the relative strong phase of the favoured and suppressed D-decay amplitudes
Bool_t clone() const
Check whether is a clone or not.
File containing declaration of LauComplex class.
LauParameter * gamma_
the relative CP-violating (weak) phase of the b -> u and b -> c amplitudes
LauParameter * createClone(Double_t constFactor=1.0)
Method to create a clone from the parent parameter using the copy constructor.
LauParameter * deltaB_
the relative CP-conserving (strong) phase of the b -> u and b -> c amplitudes
void blindParameter(const TString &blindingString, const Double_t width)
Blind the parameter.
Double_t blindingWidth() const
Obtain the Gaussian width.
Definition: LauBlind.hh:82
Double_t abs2() const
Obtain the square of the absolute value of the complex number.
Definition: LauComplex.hh:260
Bool_t fixed() const
Check whether the parameter is fixed or floated.
LauParameter()
Default constructor.
Definition: LauParameter.cc:40
LauParameter * x_
The real part of the b -> c amplitude.
File containing LauConstants namespace.
const TString & name() const
The parameter name.
const TString & blindingString() const
Obtain the blinding string.
Definition: LauBlind.hh:76
Double_t error() const
The error on the parameter.
static TRandom * getRandomiser()
Access the randomiser.
Class for defining a complex coefficient useful for extracting the CKM angle gamma from B -> D h h Da...
void updatePull()
Call to update the bias and pull values.
File containing declaration of LauPrint class and LauOutputLevel enum.
CloneOption
Options for cloning operation.
static Double_t maxRealImagPart_
Maximum allowed value of real/imaginary part parameters.
Class to define various output print commands.
Definition: LauPrint.hh:54
LauParameter * rB_
the magnitude of the ratio of the b -> u and b -> c amplitudes
static Double_t minMagnitude_
Minimum allowed value of magnitude parameters.
Double_t genValue() const
The value generated for the parameter.
static Double_t minRealImagPart_
Minimum allowed value of real/imaginary part parameters.
void printFormat(std::ostream &stream, Double_t value) const
Method to choose the printing format to a specified level of precision.
Definition: LauPrint.cc:43
virtual void finaliseValues()
Make sure values are in "standard" ranges, e.g. phases should be between -pi and pi.
LauParameter * rD_
the magnitude of the ratio of the favoured and suppressed D-decay amplitudes
void updateAmplitudes()
Update the amplitudes based on the new values of the parameters.
LauComplex nonCPPart_
The b -> c part of the complex coefficient.
virtual const LauComplex & particleCoeff()
Retrieve the complex coefficient for a particle.
static Double_t minPhase_
Minimum allowed value of phase parameters.
LauComplex antiparticleCoeff_
The antiparticle complex coefficient.
Double_t initValue() const
The initial value of the parameter.
const Double_t twoPi
Two times Pi.
virtual std::vector< LauParameter * > getParameters()
Retrieve the parameters of the coefficient, e.g. so that they can be loaded into a fit.