ROOTPWA
testAmplitude.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2010
4 //
5 // This file is part of rootpwa
6 //
7 // rootpwa is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // rootpwa is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with rootpwa. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //-------------------------------------------------------------------------
22 // File and Version Information:
23 // $Rev:: $: revision of last commit
24 // $Author:: $: author of last commit
25 // $Date:: $: date of last commit
26 //
27 // Description:
28 // basic test program for amplitude classes
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <fstream>
39 
40 #include <boost/progress.hpp>
41 
42 #include "TVector3.h"
43 #include "TLorentzRotation.h"
44 #include "TChain.h"
45 #include "TClonesArray.h"
46 #include "TStopwatch.h"
47 #include "TFile.h"
48 #include "TH1.h"
49 #include "TH2.h"
50 #include "TSystem.h"
51 
52 #ifdef USE_PWA2000
53 #include "Vec.h"
54 #include "lorentz.h"
55 #include "keyfile.h"
56 #include "event.h"
57 #include "wave.h"
58 extern particleDataTable PDGtable;
59 extern wave gWave;
60 #endif
61 
62 #include "mathUtils.hpp"
63 #include "reportingUtilsRoot.hpp"
64 #include "conversionUtils.hpp"
65 #include "particleDataTable.h"
66 #include "diffractiveDissVertex.h"
67 #include "massDependence.h"
68 #include "waveDescription.h"
69 #include "isobarAmplitude.h"
70 #include "isobarHelicityAmplitude.h"
71 #include "isobarCanonicalAmplitude.h"
72 #include "evtTreeHelper.h"
73 
74 
75 using namespace std;
76 using namespace boost;
77 using namespace rpwa;
78 
79 
80 int
81 main(int argc, char** argv)
82 {
83  printCompilerInfo();
84  printSvnVersion();
85 
86  rpwa::particleDataTable& pdt = rpwa::particleDataTable::instance();
87  pdt.readFile();
88  TStopwatch timer;
89 
90  if (0) {
91  isobarDecayVertex::setDebug(true);
92  decayTopology::setDebug(true);
93  isobarDecayTopology::setDebug(true);
94  massDependence::setDebug(true);
95  diffractiveDissVertex::setDebug(true);
96  isobarAmplitude::setDebug(true);
97  isobarHelicityAmplitude::setDebug(true);
98  isobarCanonicalAmplitude::setDebug(true);
99  waveDescription::setDebug(true);
100  }
101 
102  if (0) {
103 #ifdef USE_PWA2000
104  {
105  fourVec p(2, threeVec(0.5, 0.75, 1));
106  threeVec n = threeVec(0, 0, 1) / p.V();
107  cout << "before = " << n << " " << p << endl;
108  rotation R;
109  lorentzTransform L1;
110  L1.set(R.set(n.phi(), n.theta() - M_PI / 2, -M_PI / 2));
111  n *= R;
112  p *= L1;
113  cout << "L1 -> " << n << " " << p << endl;
114  lorentzTransform L2;
115  L2.set(R.set(0, signof(p.x()) * p.theta(), 0));
116  p *= L2;
117  cout << "L2 -> " << p << endl;
118  lorentzTransform L3;
119  L3.set(p);
120  p *= L3;
121  cout << "L3 -> " << p << endl;
122 
123  matrix<double> X(4, 4);
124  X = L3 * (L2 * L1);
125  lorentzTransform L(X);
126  p = fourVec(2, threeVec(0.5, 0.75, 1));
127  p *= L;
128  cout << "L -> " << p << endl;
129  }
130 #else
131  printWarn << "code disabled, because compilation of PWA2000 is disabled" << endl;
132 #endif
133 
134  {
135  TLorentzVector p(0.5, 0.75, 1, 2);
136  TVector3 n = TVector3(0, 0, 1).Cross(p.Vect());
137  TRotation R1;
138  R1.RotateZ(-n.Phi());
139  R1.RotateY(piHalf - n.Theta());
140  R1.RotateZ(piHalf);
141  n *= R1;
142  p *= R1;
143  cout << "R1 -> " << n << " " << p << endl;
144  // rotate about yHfAxis so that daughter momentum is along z-axis
145  TRotation R2;
146  R2.RotateY(-signum(p.X()) * p.Theta());
147  p *= R2;
148  cout << "R2 -> " << p << endl;
149  // boost into daughter RF
150  TLorentzRotation L3(-p.BoostVector());
151  cout << "L3 -> " << L3 * p << endl;
152 
153  R1.Transform(R2);
154  TLorentzRotation L(R1);
155  L.Boost(-p.BoostVector());
156  p = TLorentzVector(0.5, 0.75, 1, 2);
157  p *= L;
158  cout << "L -> " << p << endl;
159  }
160 
161  {
162  particlePtr X = createParticle("X");
163  TLorentzVector p(0.5, 0.75, 1, 2);
164  X->setLzVec(p);
165  isobarHelicityAmplitude amp;
166  TLorentzRotation L = amp.hfTransform(X->lzVec());
167  cout << "!!! L -> " << L * p << endl;
168  }
169  }
170 
171  if (0) {
172  TLorentzVector beam(1, 0.5, 180, 182);
173  TLorentzVector X (0.5, 0.75, 1, 3);
174  isobarHelicityAmplitude amp;
175  TLorentzRotation L = amp.gjTransform(beam, X);
176  cout << "!!! L -> " << L * X << endl;
177  }
178 
179  if (0) {
180  // define final state particles
181  particlePtr pi0 = createParticle("pi-", 0);
182  particlePtr pi1 = createParticle("pi+", 0);
183  particlePtr pi2 = createParticle("pi-", 1);
184  particlePtr pi3 = createParticle("pi+", 1);
185  particlePtr pi4 = createParticle("pi-", 2);
186  // define isobars
187  particlePtr sigma = createParticle("sigma");
188  particlePtr a1 = createParticle("a1(1260)+");
189  particlePtr f1 = createParticle("f1(1285)");
190  // define X-system
191  // 2I G 2J P C 2M refl
192  particlePtr X = createParticle("X-", 2, -1, 4, +1, 0, 2, +1);
193  // define production vertex
194  particlePtr beam = createParticle("pi-");
195  particlePtr target = createParticle("p+");
196  diffractiveDissVertexPtr prodVert = createDiffractiveDissVertex(beam, target, X);
197  // define vertices
198  massDependencePtr massDep = createRelativisticBreitWigner();
199  isobarDecayVertexPtr vert0 = createIsobarDecayVertex(X, pi4, f1, 2, 2);
200  isobarDecayVertexPtr vert1 = createIsobarDecayVertex(f1, pi2, a1, 2, 2, massDep);
201  isobarDecayVertexPtr vert2 = createIsobarDecayVertex(a1, pi3, sigma, 2, 0, massDep);
202  isobarDecayVertexPtr vert3 = createIsobarDecayVertex(sigma, pi0, pi1, 0, 0, massDep);
203  // set Lorentz vectors
204  beam->setLzVec(TLorentzVector(0.104385398, 0.0132061851, 189.987978, 189.988058));
205  pi0->setLzVec(TLorentzVector(-0.0761465106, -0.116917817, 5.89514709, 5.89844947));
206  pi1->setLzVec(TLorentzVector(-0.0244305532, -0.106013023, 30.6551865, 30.6556973));
207  pi2->setLzVec(TLorentzVector(0.000287952441, 0.10263611, 3.95724077, 3.96103114));
208  pi3->setLzVec(TLorentzVector(0.0299586212, 0.176440177, 115.703054, 115.703277));
209  pi4->setLzVec(TLorentzVector(0.176323963, -0.0985753246, 30.9972271, 30.9981995));
210  // build graph
211  vector<isobarDecayVertexPtr> decayVertices;
212  decayVertices.push_back(vert3);
213  decayVertices.push_back(vert1);
214  decayVertices.push_back(vert2);
215  decayVertices.push_back(vert0);
216  vector<particlePtr> fsParticles;
217  fsParticles.push_back(pi0);
218  fsParticles.push_back(pi1);
219  fsParticles.push_back(pi2);
220  fsParticles.push_back(pi3);
221  fsParticles.push_back(pi4);
222  isobarDecayTopologyPtr topo = createIsobarDecayTopology(prodVert, decayVertices, fsParticles);
223  // topo->checkTopology();
224  // topo->checkConsistency();
225  topo->fillKinematicsDataCache();
226  isobarHelicityAmplitude amp(topo);
227  cout << topo;
228  complex<double> decayAmp = amp.amplitude();
229  cout << "!!! decay amplitude = " << decayAmp << endl;
230 
231  if (1) { // compare to PWA2000
232 #ifdef USE_PWA2000
233  PDGtable.initialize("../../keyfiles/key5pi/pdgTable.txt");
234  event ev;
235  ifstream eventData("testEvents.evt");
236  ev.setIOVersion(1);
237  if (not (eventData >> ev).eof()) {
238  keyfile key;
239  key.open("testAmplitude.key");
240  complex<double> pwa2000amp;
241  key.run(ev, pwa2000amp, true);
242  key.rewind();
243  key.close();
244  cout << "!!! PWA2000 amplitude = " << pwa2000amp << endl;
245  cout << "!!! my amplitude = " << decayAmp << " vs. PWA2000 amplitude = " << pwa2000amp << ", "
246  << "delta = " << decayAmp - pwa2000amp << endl;
247  }
248 #else
249  printWarn << "code disabled, because compilation of PWA2000 is disabled" << endl;
250 #endif
251  }
252  }
253 
254  if (0) {
255  // define final state particles
256  particlePtr pi0 = createParticle("pi-", 0);
257  particlePtr pi1 = createParticle("pi+", 0);
258  // define X-system
259  // 2I G 2J P C 2M
260  particlePtr X = createParticle("X0", 2, +1, 2, -1, -1, 0);
261  // define production vertex
262  particlePtr beam = createParticle("pi-");
263  particlePtr target = createParticle("p+");
264  diffractiveDissVertexPtr prodVert = createDiffractiveDissVertex(beam, target, X);
265  // define vertices and final-state particles
266  isobarDecayVertexPtr vert0 = createIsobarDecayVertex(X, pi0, pi1, 2, 0);
267  vector<isobarDecayVertexPtr> decayVertices;
268  decayVertices.push_back(vert0);
269  vector<particlePtr> fsParticles;
270  fsParticles.push_back(pi0);
271  fsParticles.push_back(pi1);
272  isobarDecayTopologyPtr topo = createIsobarDecayTopology(prodVert, decayVertices, fsParticles);
273  topo->checkTopology();
274  topo->checkConsistency();
275  isobarAmplitude::setDebug (true);
276  isobarHelicityAmplitude::setDebug (true);
277  isobarCanonicalAmplitude::setDebug(true);
278  isobarAmplitudePtr amp[2] = {createIsobarHelicityAmplitude (topo),
279  createIsobarCanonicalAmplitude(topo)};
280  for (unsigned int i = 0; i < 2; ++i) {
281  amp[i]->enableBoseSymmetrization(false);
282  amp[i]->enableReflectivityBasis (false);
283  cout << *(amp[i]);
284  }
285  TChain chain("rootPwaEvtTree");
286  chain.Add("../../../massBins/2004/Q3PiData/template.both/1260.1300/1260.1300.root");
287  chain.GetListOfFiles()->ls();
288  vector<complex<double> > ampValues[2];
289 
290  TClonesArray* prodKinPartNames = 0;
291  TClonesArray* decayKinPartNames = 0;
292  if (not getParticleNamesFromRootFile(*(chain.GetFile()), prodKinPartNames, decayKinPartNames,
293  "rootPwaEvtTree"))
294  cout << "cannot find production and/or decay kinematics particle names "
295  << "in input file '" << chain.GetFile()->GetName() << "'." << endl;
296  for (unsigned int i = 0; i < 2; ++i)
297  processTree(chain, *prodKinPartNames, *decayKinPartNames, amp[i], ampValues[i], 2);
298  for (unsigned int i = 0; i < ampValues[0].size(); ++i)
299  cout << "amplitude[" << i << "] = " << ampValues[0][i] << " vs. " << ampValues[1][i] << "; "
300  << "ratio = " << ampValues[0][i] / ampValues[1][i] << endl;
301  }
302 
303  if (1) {
304  const long int maxNmbEvents = 10000;
305  //const long int maxNmbEvents = 1;
306 
307  // const string newKeyFileName = "test.key";
308  // const string oldKeyFileName = "testAmplitude.key";
309 
310  // const string newKeyFileName = "../../keyfiles/key3pi/SET1_new/1-0-+0+rho770_11_pi-.key";
311  // const string oldKeyFileName = "../../keyfiles/key3pi/SET1/1-0-+0+rho770_11_pi-.key";
312  // const string evtInFileName = "../../../massBins/2004/Q3PiData/template.both/1260.1300/1260.1300.evt";
313  // const string rootInFileName = "../../../massBins/2004/Q3PiData/template.both/1260.1300/1260.1300.root";
314 
315  // waves w/o isospin symmetrization
316 
317  // // rel. delta = (1.4218332033726580e-09, 1.4065747228912715e-09)
318  // // rms 1.44e-9, 8.34e-10
319  // const string newKeyFileName = "test5pi/charly/nosym/1-0-00+f01500=sigma_00_sigma_00_pi-.key";
320  // const string oldKeyFileName = "test5pi/sebastian/nosym/1-0-+0+pi-_00_f01500=sigma_0_sigma.key";
321 
322  // // rel. delta = (1.3390541006872489e-09, 1.3820297766255656e-09)
323  // // rms 1.15e-9, 1.19e-9
324  // const string newKeyFileName = "test5pi/charly/nosym/1-1+00+sigma_22_a21320-=rho770_21_pi-.key";
325  // const string oldKeyFileName = "test5pi/sebastian/nosym/1-1++0+sigma_22_a21320=pi-_2_rho770.key";
326 
327  // // rel. delta = (1.3969833147493598e-09, 1.3631710002893191e-09)
328  // // rms 1.71e-9, 1.46e-9
329  // const string newKeyFileName = "test5pi/charly/nosym/1-2-00+rho770_02_a21320-=rho770_21_pi-.key";
330  // const string oldKeyFileName = "test5pi/sebastian/nosym/1-2-+0+rho770_02_a21320=pi-_2_rho770.key";
331 
332  // rel. delta = (4.2414900268409422e-01, -3.7760164904837606e-01)
333  // rel. delta = (-1.1363857545500716e-01, 1.1934564489491958e-03)
334  // rel. delta = (-2.1179976175700090e-05, 3.5535687503020203e-06)
335  // rel. delta = (-3.4768493036075876e-09, 3.5819563781228545e-07) noBose noRefl
336  // rel. delta = (-1.3114336807010326e-07, 3.4071135025521628e-07) full sym
337  // rel. delta = (-4.2926706747269698e-05, -6.5415181989729538e-06) full sym; orig .evt
338  // rel. delta = (-7.1429252866706858e-13, -8.0913804091617701e-12) noBose noRefl; fixed PDG table
339  // rel. delta = (1.3880163165229989e-09, 1.3994113453639803e-09) full sym
340  // rel. delta = (-4.2794189442389053e-05, -6.8808316461363639e-06) full sym; orig .evt
341  // rms 1.21e-9, 2.92e-10
342  // const string newKeyFileName = "test5pi/charly/nosym/1-2-00+sigma_20_pi1800-=sigma_00_pi-.key";
343  // const string oldKeyFileName = "test5pi/sebastian/nosym/1-2-+0+sigma_20_pi1800=pi-_0_sigma.key";
344 
345  // waves with isospin symmetrization
346 
347  // rel. delta = (1.4271813795388615e-09, 1.3774171038156871e-09)
348  // rel. delta = (1.4319815006796074e-09, 1.4120124973618052e-09)
349  // rms 1.51e-9, 1.91e-9
350  // phi = 0, R = 1 ---> 1 / sqrt(2) * (a1 + a2)
351  // const string newKeyFileName = "test5pi/charly/sym/1-1+00+rho1700=a11260-=rho770_01_pi-_01_pi+_01_pi-.key";
352  // const string oldKeyFileName = "test5pi/sebastian/sym/1-1++0+pi-_01_rho1700=a11269=pi-_0_rho770_0_pi+.key";
353  const string newKeyFileName = "test5pi/charly/sym/1-1+00+rho1700=a11260+=rho770_01_pi+_01_pi-_01_pi-.key";
354  const string oldKeyFileName = "test5pi/sebastian/sym/1-1++0+pi-_01_rho1700=a11269=pi+_0_rho770_0_pi-.key";
355 
356  const string evtInFileName = "test5pi/1900.1960.genbod.regen.evt";
357  const string rootInFileName = "test5pi/1900.1960.genbod.root";
358  // const string evtInFileName = "test5pi/oneEvent.evt";
359  // const string rootInFileName = "test5pi/oneEvent.root";
360 
361  // decayTopology::setDebug(true);
362  // isobarDecayTopology::setDebug(true);
363  //massDependence::setDebug(true);
364  //isobarAmplitude::setDebug(true);
365  //isobarHelicityAmplitude::setDebug(true);
366 
367  waveDescription waveDesc;
368  isobarAmplitudePtr amp;
369  if (waveDesc.parseKeyFile(newKeyFileName) and waveDesc.constructAmplitude(amp)) {
370  isobarDecayTopologyPtr topo = amp->decayTopology();
371  printInfo << *amp;
372  amp->init();
373 
374  // read data from tree
375  const string& inTreeName = "rootPwaEvtTree";
376  const string& prodKinPartNamesObjName = "prodKinParticles";
377  const string& prodKinMomentaLeafName = "prodKinMomenta";
378  const string& decayKinPartNamesObjName = "decayKinParticles";
379  const string& decayKinMomentaLeafName = "decayKinMomenta";
380  vector<complex<double> > myAmps;
381  // open input file
382  vector<TTree*> inTrees;
383  TClonesArray* prodKinPartNames = 0;
384  TClonesArray* decayKinPartNames = 0;
385  {
386  vector<string> rootFileNames(1, rootInFileName);
387  vector<string> evtFileNames;
388  if (not openRootEvtFiles(inTrees, prodKinPartNames, decayKinPartNames,
389  rootFileNames, evtFileNames,
390  inTreeName, prodKinPartNamesObjName, prodKinMomentaLeafName,
391  decayKinPartNamesObjName, decayKinMomentaLeafName, true)) {
392  printErr << "problems opening input files. aborting." << endl;
393  exit(1);
394  }
395  }
396  const long int nmbEventsChain = inTrees[0]->GetEntries();
397 
398  // create branch pointers and leaf variables
399  TBranch* prodKinMomentaBr = 0;
400  TBranch* decayKinMomentaBr = 0;
401  TClonesArray* prodKinMomenta = 0;
402  TClonesArray* decayKinMomenta = 0;
403 
404  // connect leaf variables to tree branches
405  inTrees[0]->SetBranchAddress(prodKinMomentaLeafName.c_str(), &prodKinMomenta, &prodKinMomentaBr );
406  inTrees[0]->SetBranchAddress(decayKinMomentaLeafName.c_str(), &decayKinMomenta, &decayKinMomentaBr);
407 
408  // loop over events
409  if (not topo->initKinematicsData(*prodKinPartNames, *decayKinPartNames)) {
410  printErr << "problems initializing input data. aborting." << endl;
411  exit(1);
412  }
413  const long int nmbEvents = ((maxNmbEvents > 0) ? min(maxNmbEvents, nmbEventsChain)
414  : nmbEventsChain);
415  progress_display progressIndicator(nmbEvents);
416  timer.Reset();
417  timer.Start();
418  for (long int eventIndex = 0; eventIndex < nmbEvents; ++eventIndex) {
419  ++progressIndicator;
420 
421  if (inTrees[0]->LoadTree(eventIndex) < 0)
422  break;
423  // read only required branches
424  prodKinMomentaBr->GetEntry (eventIndex);
425  decayKinMomentaBr->GetEntry(eventIndex);
426 
427  if (not prodKinMomenta or not decayKinMomenta) {
428  printWarn << "at least one data array is null pointer: "
429  << "prodKinMomenta = " << prodKinMomenta << ", "
430  << "decayKinMomenta = " << decayKinMomenta << ". "
431  << "skipping event." << endl;
432  continue;
433  }
434 
435  if (topo->readKinematicsData(*prodKinMomenta, *decayKinMomenta)) {
436  myAmps.push_back((*amp)());
437  if ((myAmps.back().real() == 0) or (myAmps.back().imag() == 0))
438  printWarn << "event " << eventIndex << ": " << myAmps.back() << endl;
439  topo->productionVertex()->productionAmp();
440  }
441  } // event loop
442 
443  timer.Stop();
444  printSucc << "read " << myAmps.size() << " events from file(s) "
445  << "'" << rootInFileName << "' and calculated amplitudes" << endl;
446  cout << "needed ";
447  printInfo << "myAmps[0] = " << maxPrecisionDouble(myAmps[0]) << endl;
448  timer.Print();
449 
450  if (1) { // compare to PWA2000
451 #ifdef USE_PWA2000
452  PDGtable.initialize("../../keyfiles/key5pi/pdgTable.txt");
453  ifstream eventData(evtInFileName.c_str());
454  keyfile key;
455  event ev;
456  key.open(oldKeyFileName);
457  ev.setIOVersion(1);
458  timer.Reset();
459  timer.Start();
460  long int countEvent = 0;
461  vector<complex<double> > pwa2kAmps;
462  const unsigned int logInterval = 1000;
463  while ((countEvent < maxNmbEvents) and (not (eventData >> ev).eof())) {
464  complex<double> pwa2kamp;
465  key.run(ev, pwa2kamp, true);
466  pwa2kAmps.push_back(pwa2kamp);
467  key.rewind();
468  ++countEvent;
469  if (countEvent % logInterval == 0)
470  cout << " " << countEvent << endl;
471  }
472  timer.Stop();
473  printSucc << "read " << pwa2kAmps.size() << " events from file(s) "
474  << "'" << evtInFileName << "' and calculated amplitudes" << endl;
475  cout << "needed ";
476  timer.Print();
477  printInfo << "myAmps[0] = " << maxPrecisionDouble(myAmps[0]) << " vs. pwa2kAmps[0] = "
478  << maxPrecisionDouble(pwa2kAmps[0]) << ", abs. delta = "
479  << maxPrecisionDouble(myAmps[0] - pwa2kAmps[0]) << ", rel. delta = "
480  << "(" << maxPrecision((myAmps[0].real() - pwa2kAmps[0].real()) / myAmps[0].real())
481  << ", " << maxPrecision((myAmps[0].imag() - pwa2kAmps[0].imag()) / myAmps[0].imag())
482  << ")" << endl;
483 
484  if (1) {
485  const string outFileName = "testAmplitudeDiff.root";
486  printInfo << "writing comparison plots to " << outFileName << endl;
487  TFile* f = TFile::Open(outFileName.c_str(), "RECREATE");
488  TH1D* hMyAmpsReal = new TH1D("hMyAmpsReal", "hMyAmpsReal;Event Number;#Rgothic[Amplitude]", myAmps.size(), -0.5, myAmps.size() - 0.5);
489  TH1D* hMyAmpsImag = new TH1D("hMyAmpsImag", "hMyAmpsImag;Event Number;#Jgothic[Amplitude]", myAmps.size(), -0.5, myAmps.size() - 0.5);
490  TH1D* hPwa2kAmpsReal = new TH1D("hPwa2kAmpsReal", "hPwa2kAmpsReal;Event Number;#Rgothic[Amplitude]", pwa2kAmps.size(), -0.5, pwa2kAmps.size() - 0.5);
491  TH1D* hPwa2kAmpsImag = new TH1D("hPwa2kAmpsImag", "hPwa2kAmpsImag;Event Number;#Jgothic[Amplitude]", pwa2kAmps.size(), -0.5, pwa2kAmps.size() - 0.5);
492  TH1D* hDiffReal = new TH1D("hDiffReal", "hDiffReal;#Rgothic[Amplitude] Difference;Count", 100000, -1e-7, 1e-7);
493  TH1D* hDiffImag = new TH1D("hDiffImag", "hDiffImag;#Jgothic[Amplitude] Difference;Count", 100000, -1e-7, 1e-7);
494  TH2D* hCorrReal = new TH2D("hCorrReal", "hCorrReal;#Rgothic[My Amp];#Rgothic[PWA2000 Amp]", 1000, -2, 2, 1000, -2, 2);
495  TH2D* hCorrImag = new TH2D("hCorrImag", "hCorrImag;#Jgothic[My Amp];#Jgothic[PWA2000 Amp]", 1000, -2, 2, 1000, -2, 2);
496  for (unsigned int i = 0; i < myAmps.size(); ++i) {
497  hMyAmpsReal->SetBinContent (i + 1, myAmps[i].real());
498  hMyAmpsImag->SetBinContent (i + 1, myAmps[i].imag());
499  hPwa2kAmpsReal->SetBinContent(i + 1, pwa2kAmps[i].real());
500  hPwa2kAmpsImag->SetBinContent(i + 1, pwa2kAmps[i].imag());
501  hDiffReal->Fill((pwa2kAmps[i].real() - myAmps[i].real()) / pwa2kAmps[i].real());
502  hDiffImag->Fill((pwa2kAmps[i].imag() - myAmps[i].imag()) / pwa2kAmps[i].imag());
503  // hDiffReal->Fill(pwa2kAmps[i].real() - myAmps[i].real());
504  // hDiffImag->Fill(pwa2kAmps[i].imag() - myAmps[i].imag());
505  hCorrReal->Fill(myAmps[i].real(), pwa2kAmps[i].real());
506  hCorrImag->Fill(myAmps[i].imag(), pwa2kAmps[i].imag());
507  }
508  f->Write();
509  f->Close();
510  }
511 #else
512  printWarn << "code disabled, because compilation of PWA2000 is disabled" << endl;
513 #endif
514  } // compare to PWA2000
515  } // parsing of key file successful
516 
517  }
518 }
Generated by   doxygen 1.8.1.2