ROOTPWA
testIsospinSym.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 // test program for amplitude isospin symmetrization
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <fstream>
39 #include <cassert>
40 
41 #include <boost/progress.hpp>
42 
43 #include "TTree.h"
44 #include "TClonesArray.h"
45 #include "TFile.h"
46 #include "TH1.h"
47 #include "TH2.h"
48 #include "TStopwatch.h"
49 
50 #ifdef USE_PWA2000
51 #include "keyfile.h"
52 #include "event.h"
53 extern particleDataTable PDGtable;
54 #endif
55 
56 #include "mathUtils.hpp"
57 #include "fileUtils.hpp"
58 #include "particleDataTable.h"
59 #include "waveDescription.h"
60 #include "isobarHelicityAmplitude.h"
61 #include "evtTreeHelper.h"
62 
63 
64 using namespace std;
65 using namespace boost;
66 using namespace rpwa;
67 
68 
69 long int
70 calcNewAmps(const string& rootInFileName,
71  const string& keyFileName,
72  vector<complex<double> >& amps,
73  const long int maxNmbEvents)
74 {
75  waveDescription waveDesc;
76  isobarAmplitudePtr amp;
77  if (not waveDesc.parseKeyFile(keyFileName) or not waveDesc.constructAmplitude(amp)) {
78  printErr << "problems constructing amplitude. returning 0." << endl;
79  return 0;
80  }
81  isobarDecayTopologyPtr topo = amp->decayTopology();
82  printInfo << *amp;
83  amp->init();
84 
85  // open input file
86  vector<TTree*> inTrees;
87  TClonesArray* prodKinPartNames = 0;
88  TClonesArray* decayKinPartNames = 0;
89  const vector<string> rootFileNames(1, rootInFileName);
90  const vector<string> evtFileNames;
91  const string& inTreeName = "rootPwaEvtTree";
92  const string& prodKinPartNamesObjName = "prodKinParticles";
93  const string& prodKinMomentaLeafName = "prodKinMomenta";
94  const string& decayKinPartNamesObjName = "decayKinParticles";
95  const string& decayKinMomentaLeafName = "decayKinMomenta";
96  if (not openRootEvtFiles(inTrees, prodKinPartNames, decayKinPartNames,
97  rootFileNames, evtFileNames,
98  inTreeName, prodKinPartNamesObjName, prodKinMomentaLeafName,
99  decayKinPartNamesObjName, decayKinMomentaLeafName, true)) {
100  printErr << "problems opening input files. returning 0." << endl;
101  return 0;
102  }
103  const long int nmbEventsChain = inTrees[0]->GetEntries();
104 
105  // create branch pointers and leaf variables
106  TBranch* prodKinMomentaBr = 0;
107  TBranch* decayKinMomentaBr = 0;
108  TClonesArray* prodKinMomenta = 0;
109  TClonesArray* decayKinMomenta = 0;
110 
111  // connect leaf variables to tree branches
112  inTrees[0]->SetBranchAddress(prodKinMomentaLeafName.c_str(), &prodKinMomenta, &prodKinMomentaBr );
113  inTrees[0]->SetBranchAddress(decayKinMomentaLeafName.c_str(), &decayKinMomenta, &decayKinMomentaBr);
114 
115  // loop over events
116  if (not topo->initKinematicsData(*prodKinPartNames, *decayKinPartNames)) {
117  printErr << "problems initializing input data. returning 0." << endl;
118  return 0;
119  }
120  const long int nmbEvents = ((maxNmbEvents > 0) ? min(maxNmbEvents, nmbEventsChain)
121  : nmbEventsChain);
122  progress_display progressIndicator(nmbEvents);
123  for (long int eventIndex = 0; eventIndex < nmbEvents; ++eventIndex) {
124  ++progressIndicator;
125 
126  if (inTrees[0]->LoadTree(eventIndex) < 0)
127  break;
128  // read only required branches
129  prodKinMomentaBr->GetEntry (eventIndex);
130  decayKinMomentaBr->GetEntry(eventIndex);
131 
132  if (not prodKinMomenta or not decayKinMomenta) {
133  printWarn << "at least one data array is null pointer: "
134  << "prodKinMomenta = " << prodKinMomenta << ", "
135  << "decayKinMomenta = " << decayKinMomenta << ". "
136  << "skipping event." << endl;
137  continue;
138  }
139 
140  if (topo->readKinematicsData(*prodKinMomenta, *decayKinMomenta)) {
141  amps.push_back((*amp)());
142  if ((amps.back().real() == 0) or (amps.back().imag() == 0))
143  printWarn << "event " << eventIndex << ": " << amps.back() << endl;
144  // topo->productionVertex()->productionAmp();
145  }
146  }
147  return nmbEvents;
148 }
149 
150 
151 long int
152 calcPwa2kAmps(const string& evtInFileName,
153  const string& keyFileName,
154  vector<complex<double> >& amps,
155  const long int maxNmbEvents)
156 {
157 #ifdef USE_PWA2000
158  ifstream eventData(evtInFileName.c_str());
159  keyfile key;
160  event ev;
161  key.open(keyFileName);
162  ev.setIOVersion(1);
163  long int countEvent = 0;
164  const long int logInterval = 1000;
165  while ((countEvent < maxNmbEvents) and (not (eventData >> ev).eof())) {
166  complex<double> amp;
167  key.run(ev, amp, true);
168  amps.push_back(amp);
169  key.rewind();
170  ++countEvent;
171  if (countEvent % logInterval == 0)
172  cout << " " << countEvent << endl;
173  }
174  return countEvent;
175 #else
176  return 0;
177  printWarn << "code disabled, because compilation of PWA2000 is disabled" << endl;
178 #endif
179 }
180 
181 
182 struct symInfo {
183  string waveNames[2];
184  string symWaveName;
185  double phi;
186  double ratio;
187 };
188 
189 
190 vector<symInfo>
191 readDatFile(const string& fileName,
192  const bool debug = false)
193 {
194  ifstream datFile(fileName.c_str());
195  vector<symInfo> symInfos;
196  while (datFile.good()) {
197  symInfo s;
198  getline(datFile, s.waveNames[0]);
199  getline(datFile, s.waveNames[1]);
200  getline(datFile, s.symWaveName);
201  string phi;
202  getline(datFile, phi);
203  if ((phi == "pi") or (phi == "3.1416") or (phi == "1"))
204  s.phi = rpwa::pi;
205  else
206  s.phi = atof(phi.c_str());
207  string ratio;
208  getline(datFile, ratio);
209  s.ratio = atof(ratio.c_str());
210  symInfos.push_back(s);
211 
212  if (debug)
213  printDebug << "read " << symInfos.size() << ". symInfo from file '" << fileName << "'" << endl
214  << " name1 = '" << s.waveNames[0] << "'" << endl
215  << " name1 = '" << s.waveNames[1] << "'" << endl
216  << " symName = '" << s.symWaveName << "'" << endl
217  << " phi = " << maxPrecision(s.phi) << endl
218  << " ratio = " << maxPrecision(s.ratio) << endl;
219 
220  }
221  return symInfos;
222 }
223 
224 
225 int
226 main(int argc, char** argv)
227 {
228  printCompilerInfo();
229  printSvnVersion();
230 
231  const long int maxNmbEvents = 1000;
232  //const long int maxNmbEvents = 1;
233 
234  rpwa::particleDataTable& pdt = rpwa::particleDataTable::instance();
235  pdt.readFile();
236  TStopwatch timer;
237 
238  // waves with isospin symmetrization
239  // rel. delta = (1.4396529281641782e-09, 1.4018652315552070e-09)
240  // rms 4.88e-10, 1.56e-9
241  // const string newKeyFileName = "test5pi/charly/sym/1-1+00+rho1700=a11260-=rho770_01_pi-_01_pi+_01_pi-.key";
242  // const string newKeyFileName = "test5pi/charly/sym/1-1+00+rho1700=a11260+=rho770_01_pi+_01_pi-_01_pi-.key";
243  // rel. delta = (1.3916463537100299e-09, 1.3553984601344382e-09)
244  // rms 5.65e-10, 1.89e-9
245  const string newKeyFileName = "test5pi/charly/sym/1-1+00+f11285=a11260-=rho770_01_pi-_11_pi+_11_pi-.key";
246  // const string newKeyFileName = "test5pi/charly/sym/1-1+00+f11285=a11260+=rho770_01_pi+_11_pi-_11_pi-.key";
247  // rel. delta = (1.4267325099126538e-09, 1.4274825765880905e-09)
248  // rms 1.50e-9, 3.18e-9
249  // const string newKeyFileName = "test5pi/charly/sym/1-2-00+f21270=a11260-=rho770_01_pi-_11_pi+_02_pi-.key";
250  // rel. delta = (1.3479408700334261e-09, 1.3596883619015898e-09)
251  // rms 2.48e-9, 6.67e-10
252  // const string newKeyFileName = "test5pi/charly/sym/1-2-00+f21270=a11260-=rho770_01_pi-_11_pi+_22_pi-.key";
253  const string pwa2kKeyFileName[2] = {
254  // "test5pi/sebastian/sym/1-1++0+pi-_01_rho1700=a11269=pi+_0_rho770_0_pi-.key",
255  // "test5pi/sebastian/sym/1-1++0+pi-_01_rho1700=a11269=pi-_0_rho770_0_pi+.key"
256  // "test5pi/sebastian/sym/1-1++0+pi-_01_rho1700=a11269=pi+_0_rho770_0_pi-_noBose.key",
257  // "test5pi/sebastian/sym/1-1++0+pi-_01_rho1700=a11269=pi-_0_rho770_0_pi+_noBose.key"
258  "test5pi/sebastian/sym/1-1++0+pi-_11_f11285=pi-_11_a11269=pi+_0_rho770.key",
259  "test5pi/sebastian/sym/1-1++0+pi-_11_f11285=pi+_11_a11269=pi-_0_rho770.key"
260  // "test5pi/sebastian/sym/1-2-+0+pi-_02_f21270=pi-_1_a11269=pi+_0_rho770.key",
261  // "test5pi/sebastian/sym/1-2-+0+pi-_02_f21270=pi+_1_a11269=pi-_0_rho770.key"
262  // "test5pi/sebastian/sym/1-2-+0+pi-_22_f21270=pi-_1_a11269=pi+_0_rho770.key",
263  // "test5pi/sebastian/sym/1-2-+0+pi-_22_f21270=pi+_1_a11269=pi-_0_rho770.key"
264  };
265 
266  // read symmetrization info from .dat files
267  vector<string> datFiles = filesMatchingGlobPattern("test5pi/sebastian/sym/*.dat");
268  vector<symInfo> symInfos;
269  for (unsigned int i = 0; i < datFiles.size(); ++i) {
270  vector<symInfo> s = readDatFile(datFiles[i]);
271  symInfos.insert(symInfos.end(), s.begin(), s.end());
272  }
273  // get symmetrization parameters from .dat files
274  double phi = 0;
275  double ratio = 0;
276  for (unsigned int i = 0; i < symInfos.size(); ++i) {
277  const string waveName = fileNameNoExtFromPath(pwa2kKeyFileName[0]) + ".amp";
278  if ( (symInfos[i].waveNames[0] == waveName)
279  or (symInfos[i].waveNames[1] == waveName)) {
280  phi = symInfos[i].phi;
281  ratio = symInfos[i].ratio;
282  break;
283  }
284  }
285  printInfo << "using phi = " << maxPrecision(phi) << " and ratio = " << maxPrecision(ratio)
286  << " for symmetrization" << endl;
287 
288  const string evtInFileName = "test5pi/1900.1960.genbod.regen.evt";
289  const string rootInFileName = "test5pi/1900.1960.genbod.root";
290  // const string evtInFileName = "test5pi/oneEvent.evt";
291  // const string rootInFileName = "test5pi/oneEvent.root";
292 
293  //decayTopology::setDebug (true);
294  //isobarDecayTopology::setDebug (true);
295  //massDependence::setDebug (true);
296  //isobarAmplitude::setDebug (true);
297  //isobarHelicityAmplitude::setDebug(true);
298 
299  timer.Reset();
300  timer.Start();
301  vector<complex<double> > newAmps;
302  calcNewAmps(rootInFileName, newKeyFileName, newAmps, maxNmbEvents);
303  timer.Stop();
304  printSucc << "read " << newAmps.size() << " events from file(s) "
305  << "'" << rootInFileName << "' and calculated amplitudes" << endl;
306  cout << "needed ";
307  printInfo << "newAmps[0] = " << maxPrecisionDouble(newAmps[0]) << endl;
308  timer.Print();
309 
310 #ifdef USE_PWA2000
311  PDGtable.initialize("../../keyfiles/key5pi/pdgTable.txt");
312 
313  timer.Reset();
314  timer.Start();
315  vector<complex<double> > pwa2kAmps;
316  {
317  vector<complex<double> > amps[2];
318  calcPwa2kAmps(evtInFileName, pwa2kKeyFileName[0], amps[0], maxNmbEvents);
319  calcPwa2kAmps(evtInFileName, pwa2kKeyFileName[1], amps[1], maxNmbEvents);
320  const unsigned int nmbAmps = amps[0].size();
321  assert(amps[1].size() == nmbAmps);
322  pwa2kAmps.resize(nmbAmps);
323  complex<double> phase(ratio * cos(phi), ratio * sin(phi));
324  printInfo << "PWA2000 term[0] = " << maxPrecisionDouble(amps[0][0])
325  << ", term[1] = " << maxPrecisionDouble(amps[1][0])
326  << ", phase = " << phase << endl;
327  for (unsigned int i = 0; i < nmbAmps; ++i)
328  pwa2kAmps[i] = (1 / sqrt(1 + ratio)) * (amps[0][i] + phase * amps[1][i]);
329  }
330  timer.Stop();
331  printSucc << "read " << pwa2kAmps.size() << " events from file(s) "
332  << "'" << evtInFileName << "' and calculated amplitudes" << endl;
333  cout << "needed ";
334  timer.Print();
335  {
336  double relDiff[2] = {(pwa2kAmps[0].real() - newAmps[0].real()) / pwa2kAmps[0].real(),
337  (pwa2kAmps[0].imag() - newAmps[0].imag()) / pwa2kAmps[0].imag()};
338  for (unsigned int i = 0; i < 2; ++i) {
339  if (relDiff[i] < 1)
340  relDiff[i] += 2;
341  else if (relDiff[i] > 1)
342  relDiff[i] -= 2;
343  }
344  printInfo << "newAmps[0] = " << maxPrecisionDouble(newAmps[0]) << " vs. pwa2kAmps[0] = "
345  << maxPrecisionDouble(pwa2kAmps[0]) << ", abs. delta = "
346  << maxPrecisionDouble(newAmps[0] - pwa2kAmps[0]) << ", rel. delta = "
347  << "(" << maxPrecision(relDiff[0]) << ", " << maxPrecision(relDiff[1]) << ")" << endl;
348  }
349 
350  if (1) {
351  const string outFileName = "testAmplitudeDiff.root";
352  printInfo << "writing comparison plots to " << outFileName << endl;
353  TFile* f = TFile::Open(outFileName.c_str(), "RECREATE");
354  TH1D* hMyAmpsReal = new TH1D("hMyAmpsReal", "hMyAmpsReal;Event Number;#Rgothic[Amplitude]", newAmps.size(), -0.5, newAmps.size() - 0.5);
355  TH1D* hMyAmpsImag = new TH1D("hMyAmpsImag", "hMyAmpsImag;Event Number;#Jgothic[Amplitude]", newAmps.size(), -0.5, newAmps.size() - 0.5);
356  TH1D* hPwa2kAmpsReal = new TH1D("hPwa2kAmpsReal", "hPwa2kAmpsReal;Event Number;#Rgothic[Amplitude]", pwa2kAmps.size(), -0.5, pwa2kAmps.size() - 0.5);
357  TH1D* hPwa2kAmpsImag = new TH1D("hPwa2kAmpsImag", "hPwa2kAmpsImag;Event Number;#Jgothic[Amplitude]", pwa2kAmps.size(), -0.5, pwa2kAmps.size() - 0.5);
358  TH1D* hDiffReal = new TH1D("hDiffReal", "hDiffReal;#Rgothic[Amplitude] Difference;Count", 100000, -1e-7, 1e-7);
359  TH1D* hDiffImag = new TH1D("hDiffImag", "hDiffImag;#Jgothic[Amplitude] Difference;Count", 100000, -1e-7, 1e-7);
360  TH2D* hCorrReal = new TH2D("hCorrReal", "hCorrReal;#Rgothic[My Amp];#Rgothic[PWA2000 Amp]", 1000, -2, 2, 1000, -2, 2);
361  TH2D* hCorrImag = new TH2D("hCorrImag", "hCorrImag;#Jgothic[My Amp];#Jgothic[PWA2000 Amp]", 1000, -2, 2, 1000, -2, 2);
362  for (unsigned int i = 0; i < newAmps.size(); ++i) {
363  double relDiff[2] = {(pwa2kAmps[i].real() - newAmps[i].real()) / pwa2kAmps[i].real(),
364  (pwa2kAmps[i].imag() - newAmps[i].imag()) / pwa2kAmps[i].imag()};
365  for (unsigned int j = 0; j < 2; ++j) {
366  if (relDiff[j] < 1)
367  relDiff[j] += 2;
368  else if (relDiff[j] > 1)
369  relDiff[j] -= 2;
370  }
371  hMyAmpsReal->SetBinContent (i + 1, newAmps[i].real());
372  hMyAmpsImag->SetBinContent (i + 1, newAmps[i].imag());
373  hPwa2kAmpsReal->SetBinContent(i + 1, pwa2kAmps[i].real());
374  hPwa2kAmpsImag->SetBinContent(i + 1, pwa2kAmps[i].imag());
375  hDiffReal->Fill(relDiff[0]);
376  hDiffImag->Fill(relDiff[1]);
377  // hDiffReal->Fill(pwa2kAmps[i].real() - newAmps[i].real());
378  // hDiffImag->Fill(pwa2kAmps[i].imag() - newAmps[i].imag());
379  hCorrReal->Fill(newAmps[i].real(), pwa2kAmps[i].real());
380  hCorrImag->Fill(newAmps[i].imag(), pwa2kAmps[i].imag());
381  }
382  f->Write();
383  f->Close();
384  }
385 #else
386  printWarn << "code disabled, because compilation of PWA2000 is disabled" << endl;
387 #endif
388 }
Generated by   doxygen 1.8.1.2