ROOTPWA
testFitResult.C
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <vector>
4 #include <complex>
5 
6 #include "TTree.h"
7 #include "TMatrixT.h"
8 #include "TComplex.h"
9 
10 #include "reportingUtils.hpp"
11 #include "TFitResult.h"
12 #include "fitResult.h"
13 
14 
15 using namespace std;
16 using namespace rpwa;
17 
18 
19 #ifdef USE_TFITRESULT
20 
21 
22 void
23 testFitResult(TTree* oldTree,
24  TTree* newTree,
25  const bool verbose = true,
26  const string& oldBranchName = "fitResult",
27  const string& newBranchName = "fitResult_v2")
28 {
29 
30  const bool copyBin = false;
31 
32  TFitResult* oldResult = 0;
33  oldTree->SetBranchAddress(oldBranchName.c_str(), &oldResult);
34  oldTree->GetEntry(0);
35  fitResult* newResult = 0;
36  if (!copyBin) {
37  newTree->SetBranchAddress(newBranchName.c_str(), &newResult);
38  newTree->GetEntry(0);
39  } else
40  newResult = new fitResult(*oldResult);
41  const unsigned int n = newResult->nmbWaves();
42 
43  if (0)
44  cout << *newResult << endl;
45 
46  if (1) {
47  complex<double> maxDelta = 0;
48  for (unsigned int i = 0; i < n; ++i)
49  for (unsigned int j = 0; j < n; ++j) {
50  const complex<double> oldVal = oldResult->spinDensityMatrixElem(i, j);
51  const complex<double> newVal = newResult->spinDensityMatrixElem(i, j);
52  const complex<double> delta = oldVal - newVal;
53  maxDelta.real() = (fabs(maxDelta.real()) < fabs(delta.real())) ? delta.real() : maxDelta.real();
54  maxDelta.imag() = (fabs(maxDelta.imag()) < fabs(delta.imag())) ? delta.imag() : maxDelta.imag();
55  if (verbose)
56  cout << "spinDensityMatrixElem(" << newResult->waveName(i) << ", " << newResult->waveName(j) << "): "
57  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << delta << endl;
58  }
59  cout << "spinDensityMatrixElem() max. deviation = " << maxDelta << endl << endl;
60  }
61 
62  if (1) {
63  double maxDelta = 0;
64  for (unsigned int i = 0; i < n; ++i) {
65  const double oldVal = oldResult->intensity(i);
66  const double newVal = newResult->intensity(i);
67  const double delta = oldVal - newVal;
68  maxDelta = (fabs(maxDelta) < fabs(delta)) ? delta : maxDelta;
69  if (verbose)
70  cout << "intensity(" << newResult->waveName(i) << "): "
71  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << oldVal - newVal << endl;
72  }
73  cout << "intensity() max. deviation = " << maxDelta << endl << endl;
74  }
75 
76  if (1) {
77  double maxDelta = 0;
78  for (unsigned int i = 0; i < n; ++i) {
79  const double oldVal = oldResult->intensityErr(i);
80  const double newVal = newResult->intensityErr(i);
81  const double delta = oldVal - newVal;
82  maxDelta = (fabs(maxDelta) < fabs(delta)) ? delta : maxDelta;
83  if (verbose)
84  cout << "intensityErr(" << newResult->waveName(i) << "): "
85  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << oldVal - newVal << endl;
86  }
87  cout << "intensityErr() max. deviation = " << maxDelta << endl << endl;
88  }
89 
90  if (1) {
91  double maxDelta = 0;
92  const string waveNamePatterns[] = {"", // total intensity
93  "flat",
94  "0++0-",
95  "0-+0+",
96  "1++0+",
97  "2-+0+",
98  "2++0-",
99  "2++1+",
100  "1-+0-",
101  "1-+1+",
102  "3++0+",
103  "4++1+",
104  "3-+1+",
105  "3-+1-",
106  "3-+0-"};
107  for (unsigned int i = 0; i < sizeof(waveNamePatterns) / sizeof(string); ++i) {
108  const double oldVal = oldResult->intensity(waveNamePatterns[i].c_str());
109  const double newVal = newResult->intensity(waveNamePatterns[i].c_str());
110  const double delta = oldVal - newVal;
111  maxDelta = (fabs(maxDelta) < fabs(delta)) ? delta : maxDelta;
112  if (verbose)
113  cout << "intensity(" << waveNamePatterns[i] << "): "
114  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << oldVal - newVal << endl;
115  }
116  cout << "intensity() max. deviation = " << maxDelta << endl << endl;
117  }
118 
119  if (1) {
120  double maxDelta = 0;
121  const string waveNamePatterns[] = {"", // total intensity
122  "flat",
123  "0++0-",
124  "0-+0+",
125  "1++0+",
126  "2-+0+",
127  "2++0-",
128  "2++1+",
129  "1-+0-",
130  "1-+1+",
131  "3++0+",
132  "4++1+",
133  "3-+1+",
134  "3-+1-",
135  "3-+0-"};
136  for (unsigned int i = 0; i < sizeof(waveNamePatterns) / sizeof(string); ++i) {
137  const double oldVal = oldResult->intensityErr(waveNamePatterns[i].c_str());
138  const double newVal = newResult->intensityErr(waveNamePatterns[i].c_str());
139  const double delta = oldVal - newVal;
140  maxDelta = (fabs(maxDelta) < fabs(delta)) ? delta : maxDelta;
141  if (verbose)
142  cout << "intensityErr(" << waveNamePatterns[i] << "): "
143  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << oldVal - newVal << endl;
144  }
145  cout << "intensityErr() max. deviation = " << maxDelta << endl << endl;
146  }
147 
148  if (1) {
149  double maxDelta = 0;
150  for (unsigned int i = 0; i < n; ++i)
151  for (unsigned int j = 0; j < n; ++j) {
152  const double oldVal = oldResult->phase(i, j);
153  const double newVal = newResult->phase(i, j);
154  const double delta = oldVal - newVal;
155  maxDelta = (fabs(maxDelta) < fabs(delta)) ? delta : maxDelta;
156  if (verbose)
157  cout << "phase(" << newResult->waveName(i) << ", " << newResult->waveName(j) << "): "
158  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << delta << endl;
159  }
160  cout << "phase() max. deviation = " << maxDelta << endl << endl;
161  }
162 
163  if (1) {
164  double maxDelta = 0;
165  for (unsigned int i = 0; i < n; ++i)
166  for (unsigned int j = 0; j < n; ++j) {
167  const double oldVal = oldResult->phaseErr(i, j);
168  const double newVal = newResult->phaseErr(i, j);
169  const double delta = oldVal - newVal;
170  maxDelta = (fabs(maxDelta) < fabs(delta)) ? delta : maxDelta;
171  if (verbose)
172  cout << "phaseErr(" << newResult->waveName(i) << ", " << newResult->waveName(j) << "): "
173  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << delta << endl;
174  }
175  cout << "phaseErr() max. deviation = " << maxDelta << endl << endl;
176  }
177 
178  if (1) {
179  double maxDelta = 0;
180  for (unsigned int i = 0; i < n; ++i)
181  for (unsigned int j = 0; j < n; ++j) {
182  const double oldVal = oldResult->coherence(i, j);
183  const double newVal = newResult->coherence(i, j);
184  const double delta = oldVal - newVal;
185  maxDelta = (fabs(maxDelta) < fabs(delta)) ? delta : maxDelta;
186  if (verbose)
187  cout << "coherence(" << newResult->waveName(i) << ", " << newResult->waveName(j) << "): "
188  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << delta << endl;
189  }
190  cout << "coherence() max. deviation = " << maxDelta << endl << endl;
191  }
192 
193  if (1) {
194  double maxDelta = 0;
195  for (unsigned int i = 0; i < n; ++i)
196  for (unsigned int j = 0; j < n; ++j) {
197  const double oldVal = oldResult->coherenceErr(i, j);
198  const double newVal = newResult->coherenceErr(i, j);
199  const double delta = oldVal - newVal;
200  maxDelta = (fabs(maxDelta) < fabs(delta)) ? delta : maxDelta;
201  if (verbose)
202  cout << "coherenceErr(" << newResult->waveName(i) << ", " << newResult->waveName(j) << "): "
203  << setprecision(12) << newVal << " vs. " << oldVal << ", delta = " << delta << endl;
204  }
205  cout << "coherenceErr() max. deviation = " << maxDelta << endl << endl;
206  }
207 }
208 
209 
210 #endif // USE_TFITRESULT