ROOTPWA
checkKeyFile.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 // checks validity of amplitude specified by key file by
29 // performing a number of consistency checks
30 //
31 //
32 // Author List:
33 // Boris Grube TUM (original author)
34 //
35 //
36 //-------------------------------------------------------------------------
37 
38 
39 #include <iostream>
40 #include <fstream>
41 #include <unistd.h>
42 #include <vector>
43 #include <complex>
44 
45 #include "TChain.h"
46 #include "TTree.h"
47 #include "TBranch.h"
48 #include "TClonesArray.h"
49 #include "TStopwatch.h"
50 
51 #include "fileUtils.hpp"
52 #include "particleDataTable.h"
53 #include "evtTreeHelper.h"
54 #include "waveDescription.h"
55 
56 
57 using namespace std;
58 using namespace rpwa;
59 
60 
61 void
62 usage(const string& progName,
63  const int errCode = 0)
64 {
65  cerr << "verify that given .key files do not violate physics" << endl
66  << endl
67  << "usage:" << endl
68  << progName
69  << " [-d test data -n max. # of events -p PDG file -t tree name -e max. diff. -l leaf names -v -h] key file(s)" << endl
70  << " where:" << endl
71  << " -d file path to file with test data (.evt or .root format)" << endl
72  << " -n # maximum number of events to read (default: all)" << endl
73  << " -p file path to particle data table file (default: ./particleDataTable.txt)" << endl
74  << " -t name name of tree in ROOT data files (default: rootPwaEvtTree)" << endl
75  << " -e # maximum deviation of amplitude ratios from 1 (default: 1E-6)" << endl
76  << " -l names semicolon separated object/leaf names in input data (default: 'prodKinParticles;prodKinMomenta;decayKinParticles;decayKinMomenta')" << endl
77  << " -v verbose; print debug output (default: false)" << endl
78  << " -h print help" << endl
79  << endl;
80  exit(errCode);
81 }
82 
83 
84 bool testAmplitude(TTree* inTree,
85  const string& keyFileName,
86  vector<string>& keyFileErrors,
87  const long int maxNmbEvents = -1,
88  const bool debug = false,
89  const double maxDelta = 1e-6,
90  const TClonesArray* prodKinPartNames = 0,
91  const TClonesArray* decayKinPartNames = 0,
92  const string& prodKinMomentaLeafName = "prodKinMomenta",
93  const string& decayKinMomentaLeafName = "decayKinMomenta")
94 {
95  // parse key file and create decay topology and amplitude instances
96  waveDescription waveDesc;
97  isobarDecayTopologyPtr decayTopo;
98  if ( not waveDesc.parseKeyFile(keyFileName)
99  or not waveDesc.constructDecayTopology(decayTopo)) {
100  printWarn << "problems constructing decay topology from key file '" << keyFileName << "'. "
101  << "skipping." << endl;
102  keyFileErrors.push_back("parsing errors");
103  return false;
104  }
105 
106  printInfo << "checking decay topology" << endl;
107  bool success = true;
108  if (not decayTopo->checkTopology()) {
109  keyFileErrors.push_back("problematic topology");
110  success = false;
111  }
112  if (not decayTopo->checkConsistency()) {
113  keyFileErrors.push_back("inconsistent decay");
114  success = false;
115  }
116  if (not success)
117  return false;
118 
119  if (not inTree or not prodKinPartNames or not decayKinPartNames)
120  return true; // no data; cannot check symmetry properties of amplitude
121 
122  printInfo << "calculating amplitudes for ";
123  if (maxNmbEvents < 0)
124  cout << "all";
125  else
126  cout << maxNmbEvents;
127  cout << " events" << endl;
128 
129  // construct amplitude
130  isobarAmplitudePtr amplitude;
131  waveDesc.constructAmplitude(amplitude, decayTopo);
132 
133 
134  // read data from tree and calculate amplitudes
135  vector<complex<double> > ampValues;
136  if (not processTree(*inTree, *prodKinPartNames, *decayKinPartNames,
137  amplitude, ampValues, maxNmbEvents,
138  prodKinMomentaLeafName, decayKinMomentaLeafName, false)) {
139  printWarn << "problems reading tree" << endl;
140  return false;
141  }
142  if (ampValues.size() < 1) {
143  printWarn << "no amplitude were calculated" << endl;
144  return false;
145  }
146 
147  // calculate amplitudes for parity transformed decay daughters
148  vector<complex<double> > ampSpaceInvValues;
149  amplitude->enableSpaceInversion(true);
150  if (not processTree(*inTree, *prodKinPartNames, *decayKinPartNames,
151  amplitude, ampSpaceInvValues, maxNmbEvents,
152  prodKinMomentaLeafName, decayKinMomentaLeafName, false)) {
153  printWarn << "problems reading tree" << endl;
154  return false;
155  }
156 
157  // calculate amplitudes for decay daughters reflected through production plane
158  vector<complex<double> > ampReflValues;
159  amplitude->enableSpaceInversion(false);
160  amplitude->enableReflection (true);
161  if (not processTree(*inTree, *prodKinPartNames, *decayKinPartNames,
162  amplitude, ampReflValues, maxNmbEvents,
163  prodKinMomentaLeafName, decayKinMomentaLeafName, false)) {
164  printWarn << "problems reading tree" << endl;
165  return false;
166  }
167 
168  if ( (ampValues.size() != ampSpaceInvValues.size())
169  or (ampValues.size() != ampReflValues.size ())) {
170  printWarn << "different number of amplitudes for space inverted "
171  << "(" << ampSpaceInvValues.size() << "), reflected "
172  << "(" << ampReflValues.size() << "), and unmodified data "
173  << "(" << ampValues.size() << ")." << endl;
174  return false;
175  }
176 
177  printInfo << "checking symmetry properties of amplitudes" << endl;
178  unsigned int countAmpZero = 0;
179  unsigned int countAmpRatioNotOk = 0;
180  unsigned int countSpaceInvEigenValNotOk = 0;
181  unsigned int countReflEigenValNotOk = 0;
182  for (unsigned int i = 0; i < ampValues.size(); ++i) {
183  // check that amplitude is non-zero
184  bool ampZero = false;
185  if (ampValues[i] == complex<double>(0, 0)) {
186  ampZero = true;
187  ++countAmpZero;
188  }
189  const unsigned int nmbDigits = numeric_limits<double>::digits10 + 1;
190  ostringstream s;
191  s.precision(nmbDigits);
192  s.setf(ios_base::scientific, ios_base::floatfield);
193  if (debug) {
194  printDebug << "amplitude " << i << ": " << endl;
195  s << " ampl. = " << ampValues[i];
196  if (ampZero)
197  s << " <! zero amplitude";
198  s << endl
199  << " ampl. space inv. = " << ampSpaceInvValues[i] << endl
200  << " ampl. refl. = " << ampReflValues [i] << endl;
201  cout << s.str();
202  }
203 
204  // check space inversion symmetry
205  const complex<double> spaceInvRatio
206  = complex<double>((ampSpaceInvValues[i].real() != 0) ?
207  ampValues[i].real() / ampSpaceInvValues[i].real() : 0,
208  (ampSpaceInvValues[i].imag() != 0) ?
209  ampValues[i].imag() / ampSpaceInvValues[i].imag() : 0);
210  bool spaceInvAmpRatioOk = true;
211  if ( (fabs(spaceInvRatio.real()) - 1 > maxDelta)
212  or (fabs(spaceInvRatio.imag()) - 1 > maxDelta)) {
213  spaceInvAmpRatioOk = false;
214  ++countAmpRatioNotOk;
215  }
216  const int spaceInvEigenValue = decayTopo->spaceInvEigenValue();
217  bool spaceInvEigenValOk = true;
218  if ( ( (spaceInvRatio.real() != 0)
219  and (spaceInvRatio.real() / fabs(spaceInvRatio.real()) != spaceInvEigenValue))
220  or ( (spaceInvRatio.imag() != 0)
221  and (spaceInvRatio.imag() / fabs(spaceInvRatio.imag()) != spaceInvEigenValue))) {
222  spaceInvEigenValOk = false;
223  ++countSpaceInvEigenValNotOk;
224  }
225  if (debug) {
226  s.str("");
227  s << "Re[ampl.] / Re[ampl.] space inv. - 1 = " << setw(23) << spaceInvRatio.real() - 1 << ", "
228  << "Im[ampl.] / Im[ampl.] space inv. - 1 = " << setw(23) << spaceInvRatio.imag() - 1;
229  cout << " " << s.str();
230  if (not spaceInvAmpRatioOk)
231  cout << " <! larger than " << maxDelta;
232  if (not spaceInvEigenValOk)
233  cout << " <! eigenvalue != " << spaceInvEigenValue;
234  cout << endl;
235  }
236 
237  // check reflection symmetry through production plane
238  const complex<double> reflRatio
239  = complex<double>((ampReflValues[i].real() != 0) ?
240  ampValues[i].real() / ampReflValues[i].real() : 0,
241  (ampReflValues[i].imag() != 0) ?
242  ampValues[i].imag() / ampReflValues[i].imag() : 0);
243  bool reflAmpRatioOk = true;
244  if ( (fabs(reflRatio.real()) - 1 > maxDelta)
245  or (fabs(reflRatio.imag()) - 1 > maxDelta)) {
246  reflAmpRatioOk = false;
247  ++countAmpRatioNotOk;
248  }
249  const int reflEigenValue = decayTopo->reflectionEigenValue();
250  bool reflEigenValOk = true;
251  if ( ( (reflRatio.real() != 0)
252  and (reflRatio.real() / fabs(reflRatio.real()) != reflEigenValue))
253  or ( (reflRatio.imag() != 0)
254  and (reflRatio.imag() / fabs(reflRatio.imag()) != reflEigenValue))) {
255  reflEigenValOk = false;
256  ++countReflEigenValNotOk;
257  }
258  if (debug) {
259  s.str("");
260  s << "Re[ampl.] / Re[ampl.] refl. - 1 = " << setw(23) << reflRatio.real() - 1 << ", "
261  << "Im[ampl.] / Im[ampl.] refl. - 1 = " << setw(23) << reflRatio.imag() - 1;
262  cout << " " << s.str();
263  if (not reflAmpRatioOk)
264  cout << " <! larger than " << maxDelta;
265  if (not reflEigenValOk)
266  cout << " <! eigenvalue != " << reflEigenValue;
267  cout << endl;
268  }
269  }
270 
271  if (countAmpZero > 0) {
272  keyFileErrors.push_back("zero amplitude");
273  success = false;
274  }
275  if (countAmpRatioNotOk > 0) {
276  stringstream s;
277  s << "amplitude deviation in symmetry check larger than " << maxDelta;
278  keyFileErrors.push_back(s.str());
279  success = false;
280  }
281  if (countSpaceInvEigenValNotOk > 0) {
282  keyFileErrors.push_back("wrong space inversion eigenvalue");
283  success = false;
284  }
285  if (countReflEigenValNotOk > 0) {
286  keyFileErrors.push_back("wrong eigenvalue for reflection through production plane");
287  success = false;
288  }
289  return success;
290 }
291 
292 
293 int
294 main(int argc,
295  char** argv)
296 {
297  printCompilerInfo();
298  printLibraryInfo ();
299  printSvnVersion ();
300  cout << endl;
301 
302  // parse command line options
303  const string progName = argv[0];
304  string dataFileName = "";
305  long int maxNmbEvents = -1;
306  string pdgFileName = "./particleDataTable.txt";
307  string inTreeName = "rootPwaEvtTree";
308  double maxDelta = 1e-6;
309  string leafNames = "prodKinParticles;prodKinMomenta;decayKinParticles;decayKinMomenta";
310  bool debug = false;
311  extern char* optarg;
312  extern int optind;
313  int c;
314  while ((c = getopt(argc, argv, "d:n:p:t:e:l:r:vh")) != -1)
315  switch (c) {
316  case 'd':
317  dataFileName = optarg;
318  break;
319  case 'n':
320  maxNmbEvents = atol(optarg);
321  break;
322  case 'p':
323  pdgFileName = optarg;
324  break;
325  case 't':
326  inTreeName = optarg;
327  break;
328  case 'e':
329  maxDelta = atof(optarg);
330  break;
331  case 'l':
332  leafNames = optarg;
333  break;
334  case 'v':
335  debug = true;
336  break;
337  case 'h':
338  default:
339  usage(progName);
340  }
341 
342  // set debug options
343  if (debug) {
344  waveDescription::setDebug(true);
345  //particleProperties::setDebug(true);
346  //isobarHelicityAmplitude::setDebug(true);
347  //massDependence::setDebug(true);
348  }
349 
350  // get key file names
351  if (optind >= argc) {
352  printErr << "you need to specify at least one key file to process. aborting." << endl;;
353  usage(progName, 1);
354  }
355  vector<string> keyFileNames;
356  while (optind < argc) {
357  const string fileName = argv[optind++];
358  keyFileNames.push_back(fileName);
359  }
360 
361  // get object and leaf names for event data
362  string prodKinPartNamesObjName, prodKinMomentaLeafName;
363  string decayKinPartNamesObjName, decayKinMomentaLeafName;
364  parseLeafAndObjNames(leafNames, prodKinPartNamesObjName, prodKinMomentaLeafName,
365  decayKinPartNamesObjName, decayKinMomentaLeafName);
366 
367  // open input file
368  TTree* inTree = 0;
369  TClonesArray* prodKinPartNames = 0;
370  TClonesArray* decayKinPartNames = 0;
371  if (dataFileName == "")
372  printInfo << "No test data file was specified (option -d). "
373  << "skipping symmetry checks of amplitude." << endl;
374  else {
375  vector<string> rootFileNames;
376  vector<string> evtFileNames;
377  const string fileExt = extensionFromPath(dataFileName);
378  if (fileExt == "root")
379  rootFileNames.push_back(dataFileName);
380  else if (fileExt == "evt")
381  evtFileNames.push_back(dataFileName);
382  else {
383  printErr << "specified input files is neither a .root nor a .evt file. aborting.";
384  usage(progName, 1);
385  }
386  vector<TTree*> inTrees;
387  if (not openRootEvtFiles(inTrees, prodKinPartNames, decayKinPartNames,
388  rootFileNames, evtFileNames,
389  inTreeName, prodKinPartNamesObjName, prodKinMomentaLeafName,
390  decayKinPartNamesObjName, decayKinMomentaLeafName, debug)) {
391  printErr << "problems opening input file(s). aborting." << endl;
392  exit(1);
393  }
394  inTree = inTrees[0];
395  }
396 
397  // initialize particle data table
398  particleDataTable::readFile(pdgFileName);
399 
400  // loop over key files
401  map<string, vector<string> > keyFileErrors; // maps error description to key files
402  unsigned int countKeyFileErr = 0;
403  for (unsigned int i = 0; i < keyFileNames.size(); ++i) {
404  cout << endl;
405  printInfo << "checking key file '" << keyFileNames[i] << "'" << endl;
406  vector<string> errors;
407  if (not testAmplitude(inTree, keyFileNames[i], errors, maxNmbEvents, debug, maxDelta,
408  prodKinPartNames, decayKinPartNames,
409  prodKinMomentaLeafName, decayKinMomentaLeafName)) {
410  ++countKeyFileErr;
411  printWarn << "key file '" << keyFileNames[i] << "' did not pass all tests. "
412  << "see summary below." << endl;
413  } else
414  printSucc << "key file '" << keyFileNames[i] << "' passed all tests" << endl;
415  // collect errors
416  for (unsigned int j = 0; j < errors.size(); ++j)
417  keyFileErrors[errors[j]].push_back(keyFileNames[i]);
418  }
419 
420  // clean up
421  if (inTree)
422  delete inTree;
423 
424  // report final result and set exit status accordingly
425  cout << endl;
426  if (countKeyFileErr == 0) {
427  printSucc << "all " << keyFileNames.size() << " keyfile(s) passed all tests" << endl;
428  return 0;
429  }
430 
431  printInfo << keyFileNames.size() - countKeyFileErr << " of " << keyFileNames.size()
432  << " keyfile(s) passed all tests" << endl;
433  if (keyFileErrors.size() > 0) {
434  printInfo << countKeyFileErr << " problematic keyfile(s):" << endl;
435  for (map<string, vector<string> >::const_iterator entry = keyFileErrors.begin();
436  entry != keyFileErrors.end(); ++entry) {
437  cout << " " << entry->first << ":" << endl;
438  for (unsigned int i = 0; i < entry->second.size(); ++i)
439  cout << " " << entry->second[i] << endl;
440  }
441  }
442  exit(1);
443 
444 }
Generated by   doxygen 1.8.1.2