MINT2
binflipChi2.cpp
Go to the documentation of this file.
1 #include "Mint/binflipChi2.h"
2 
3 using namespace MINT;
4 using namespace std;
5 
6 binflipChi2::binflipChi2(vector<complex<double> > X, vector<double> r, vector<double> tAv, vector<double> tSqAv, TH2F pHistD0,
7  TH2F pHistD0bar, TH2F nHistD0, TH2F nHistD0bar, double ReZcp, double ImZcp, double ReDz, double ImDz, double stepSize,
8  int fakeData, vector<double> Fm, vector<double> Fp, int verbosity):
9 
11  m_X(X),
12  m_r(r),
13  m_tAv(tAv),
14  m_tSqAv(tSqAv),
15  m_pHistD0(pHistD0),
16  m_pHistD0bar(pHistD0bar),
17  m_nHistD0(nHistD0),
18  m_nHistD0bar(nHistD0bar),
19  m_ReZcp("ReZcp", FitParameter::FIT, ReZcp, stepSize, 0, 0, getParSet()),
20  m_ImZcp("ImZcp", FitParameter::FIT, ImZcp, stepSize, 0, 0, getParSet()),
21  m_ReDz("ReDz", FitParameter::FIT, ReDz, stepSize, 0, 0, getParSet()),
22  m_ImDz("ImDz", FitParameter::FIT, ImDz, stepSize, 0, 0, getParSet()),
23  m_fakeData(fakeData),
24  m_Fm(Fm),
25  m_Fp(Fp),
26  m_verbosity(verbosity)
27 {
28  m_nbinsPhase = m_r.size();
29  m_nbinsTime = m_tAv.size();
30 
31  if (( m_fakeData > 0 ) && ( Fm.size() == m_nbinsPhase ) && ( Fp.size() == m_nbinsPhase )){
32  genFakeData();
33  }
34 }
35 
36 
38  delete getParSet();
39 }
40 
41 
43 
44  double chi2 = 0.;
45  int count = 0;
46  vector<vector<TGraph> > fits(2);
47 
48  fits = getFits();
49 
50  for(int b = 1; b <= m_nbinsPhase; b++){
51 
52  double* Rvals_pl = fits[0][b-1].GetY();
53  double* Rvals_mi = fits[1][b-1].GetY();
54 
55  for(int j = 1; j <= m_nbinsTime; j++){
56 
57  double R_pl = Rvals_pl[j-1];
58  double R_mi = Rvals_mi[j-1];
59  double D0_term = 0, D0bar_term = 0;
60  double D0_num = 0, D0_den = 0, D0bar_num = 0, D0bar_den = 0;
61 
62  if(( m_nHistD0.GetBinContent(j,b) != 0 )&&( m_pHistD0.GetBinContent(j,b) != 0 )){
63  D0_num = pow(m_nHistD0.GetBinContent(j,b) - m_pHistD0.GetBinContent(j,b) * R_pl, 2);
64  D0_den = pow(m_nHistD0.GetBinError(j,b), 2) + pow(m_pHistD0.GetBinError(j,b) * R_pl, 2);
65  if(D0_den != 0){
66  D0_term = D0_num/D0_den;
67  }
68  }
69 
70  if(( m_nHistD0bar.GetBinContent(j,b) != 0 )&&( m_pHistD0bar.GetBinContent(j,b) != 0 )){
71  D0bar_num = pow(m_nHistD0bar.GetBinContent(j,b) - m_pHistD0bar.GetBinContent(j,b) * R_mi, 2);
72  D0bar_den = pow(m_nHistD0bar.GetBinError(j,b), 2) + pow(m_pHistD0bar.GetBinError(j,b) * R_mi, 2);
73  if(D0bar_den != 0){
74  D0bar_term = D0bar_num/D0bar_den;
75  }
76  }
77 
78  chi2 += D0_term + D0bar_term;
79 
80  float tol = 0.01;
81  if(D0_term > tol){
82  count += 1;
83  }
84  if(D0bar_term > tol){
85  count += 1;
86  }
87  }
88  }
89 
90  if( m_verbosity > 0 ){
91  cout << "Number of terms contributing to chi2 is: " << count << endl;
92  }
93 
94  return chi2;
95 }
96 
97 
98 vector<vector<TGraph> > binflipChi2::getFits(){
99 
100  vector<vector<TGraph> > fits(2);
101 
102  complex<double> zcp(m_ReZcp.getCurrentFitVal(), m_ImZcp.getCurrentFitVal());
103  complex<double> deltaz(m_ReDz.getCurrentFitVal(), m_ImDz.getCurrentFitVal());
104 
105  for(int i : {0,1}){
106 
107  for(int b = 1; b <= m_nbinsPhase; b++){
108  TGraph temp = TGraph(m_nbinsTime);
109  fits[i].push_back( temp );
110 
111  for(int j = 1; j <= m_nbinsTime; j++){
112  double numerator = 0., denominator = 0.;
113 
114  numerator = m_r[b-1] * ( 1 + 0.25 * m_tSqAv[j-1] * ( pow(zcp,2) - pow(deltaz,2) ).real() );
115  numerator += 0.25 * m_tSqAv[j-1] * pow(abs(zcp + (double)pow(-1, i) * deltaz), 2);
116  numerator += sqrt(m_r[b-1]) * m_tAv[j-1] * ( conj(m_X[b-1]) * (zcp + (double)pow(-1, i) * deltaz) ).real();
117 
118  denominator = 1 + 0.25 * m_tSqAv[j-1] * ( pow(zcp,2) - pow(deltaz,2) ).real();
119  denominator += m_r[b-1] * 0.25 * m_tSqAv[j-1] * pow(abs(zcp + (double)pow(-1, i) * deltaz), 2);
120  denominator += sqrt(m_r[b-1]) * m_tAv[j-1] * ( m_X[b-1] * (zcp + (double)pow(-1, i) * deltaz) ).real();
121 
122  double Rval = numerator/denominator;
123  fits[i][b-1].SetPoint(j-1, m_tAv[j-1], Rval);
124 
125  }
126  }
127  }
128 
129  return fits;
130 }
131 
132 
134  TRandom3 rndm = TRandom3(m_fakeData);
135  complex<double> zcp(m_ReZcp.getCurrentFitVal(), m_ImZcp.getCurrentFitVal());
136  complex<double> deltaz(m_ReDz.getCurrentFitVal(), m_ImDz.getCurrentFitVal());
137 
138  complex<double> qoverp = sqrt( (zcp + deltaz) / (zcp - deltaz) );
139  complex<double> z = sqrt( pow(zcp,2) - pow(deltaz,2) );
140 
141  for(int i : {0,1}){
142 
143  complex<double> pqterm;
144  if ( i == 0 ){
145  pqterm = ((double)-1)*qoverp;
146  }
147  else{
148  pqterm = ((double)-1)*pow(qoverp, -1);
149  }
150 
151  for(int b = 1; b <= m_nbinsPhase; b++){
152 
153  for(int j = 1; j <= m_nbinsTime; j++){
154 
155  double pval = 0., nval = 0., rNum = 0., err = 0.;
156 
157  pval = m_Fp[b-1] * ( 1 + 0.25 * m_tSqAv[j-1] * ( pow(z, 2) ).real() );
158  pval += 0.25 * m_tSqAv[j-1] * pow(abs(z), 2) * pow(abs(pqterm), 2) * m_Fm[b-1];
159  pval += m_tAv[j-1] * sqrt(m_Fm[b-1] * m_Fp[b-1]) * (pqterm * m_X[b-1] * z).real();
160 
161  nval = m_Fm[b-1] * ( 1 + 0.25 * m_tSqAv[j-1] * ( pow(z, 2) ).real() );
162  nval += 0.25 * m_tSqAv[j-1] * pow(abs(z), 2) * pow(abs(pqterm), 2) * m_Fp[b-1];
163  nval += m_tAv[j-1] * sqrt(m_Fm[b-1] * m_Fp[b-1]) * (pqterm * conj(m_X[b-1]) * z).real();
164 
165  if( i == 0 ){
166  if( (m_pHistD0.GetBinContent(j,b) != 0) && (m_pHistD0.GetBinError(j, b) != 0) ){
167  err = pval * m_pHistD0.GetBinError(j, b) / m_pHistD0.GetBinContent(j,b);
168  }
169  else{
170  err = sqrt(pval);
171  }
172  if (m_pHistD0.GetBinContent(j,b) != 0){
173  rNum = rndm.Gaus(0, err);
174  m_pHistD0.SetBinContent(j, b, pval + rNum);
175  m_pHistD0.SetBinError(j, b, err);
176  }
177 
178  if( (m_nHistD0.GetBinContent(j,b) != 0) && (m_nHistD0.GetBinError(j, b) != 0) ){
179  err = nval * m_nHistD0.GetBinError(j, b) / m_nHistD0.GetBinContent(j,b);
180  }
181  else{
182  err = sqrt(nval);
183  }
184  if (m_nHistD0.GetBinContent(j,b) != 0){
185  rNum = rndm.Gaus(0, err);
186  m_nHistD0.SetBinContent(j, b, nval + rNum);
187  m_nHistD0.SetBinError(j, b, err);
188  }
189  }
190  else{
191  if( (m_pHistD0bar.GetBinContent(j,b) != 0) && (m_pHistD0bar.GetBinError(j, b) != 0) ){
192  err = pval * m_pHistD0bar.GetBinError(j, b) / m_pHistD0bar.GetBinContent(j,b);
193  }
194  else{
195  err = sqrt(pval);
196  }
197  if( m_pHistD0bar.GetBinContent(j,b) != 0 ){
198  rNum = rndm.Gaus(0, err);
199  m_pHistD0bar.SetBinContent(j, b, pval + rNum);
200  m_pHistD0bar.SetBinError(j, b, err);
201  }
202 
203  if( (m_nHistD0bar.GetBinContent(j,b) != 0) && (m_nHistD0bar.GetBinError(j, b) != 0) ){
204  err = nval * m_nHistD0bar.GetBinError(j, b) / m_nHistD0bar.GetBinContent(j,b);
205  }
206  else{
207  err = sqrt(nval);
208  }
209  if( m_nHistD0bar.GetBinContent(j,b) != 0 ){
210  rNum = rndm.Gaus(0, err);
211  m_nHistD0bar.SetBinContent(j, b, nval + rNum);
212  m_nHistD0bar.SetBinError(j, b, err);
213  }
214  }
215  }
216  }
217  }
218 }
219 
220 
221 TGraph binflipChi2::getFit(int i, int b){
222 
223  TGraph fit = TGraph(m_nbinsTime);
224  fit = getFits()[i][b];
225  return fit;
226 
227 }
FitParameter m_ReZcp
Definition: binflipChi2.h:27
binflipChi2(vector< complex< double > > X, vector< double > r, vector< double > tAv, vector< double > tSqAv, TH2F pHistD0, TH2F pHistD0bar, TH2F nHistD0, TH2F nHistD0bar, double ReZcp, double ImZcp, double ReDz, double ImDz, double stepSize, int fakeData=0, vector< double > Fm=vector< double >(), vector< double > Fp=vector< double >(), int verbosity=0)
Definition: binflipChi2.cpp:6
double getVal()
Definition: binflipChi2.cpp:42
vector< double > m_tSqAv
Definition: binflipChi2.h:22
vector< double > m_Fm
Definition: binflipChi2.h:32
TGraph getFit(int i, int b)
void genFakeData()
TH2F m_nHistD0
Definition: binflipChi2.h:25
TH2F m_pHistD0
Definition: binflipChi2.h:23
FitParameter m_ImDz
Definition: binflipChi2.h:30
TH2F m_nHistD0bar
Definition: binflipChi2.h:26
int m_nbinsPhase
Definition: binflipChi2.h:17
vector< vector< TGraph > > getFits()
Definition: binflipChi2.cpp:98
MinuitParameterSet * getParSet()
Definition: Minimisable.cpp:25
vector< double > m_Fp
Definition: binflipChi2.h:33
FitParameter m_ReDz
Definition: binflipChi2.h:29
int m_fakeData
Definition: binflipChi2.h:31
virtual double getCurrentFitVal() const
int m_verbosity
Definition: binflipChi2.h:34
vector< double > m_r
Definition: binflipChi2.h:20
int m_nbinsTime
Definition: binflipChi2.h:18
vector< complex< double > > m_X
Definition: binflipChi2.h:19
FitParameter m_ImZcp
Definition: binflipChi2.h:28
TH2F m_pHistD0bar
Definition: binflipChi2.h:24
vector< double > m_tAv
Definition: binflipChi2.h:21