ROOTPWA
ampDiff.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 // compares two amplitude files and calculates differences
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <cstdlib>
39 #include <iostream>
40 #include <fstream>
41 #include <unistd.h>
42 #include <vector>
43 #include <complex>
44 
45 #include "TFile.h"
46 #include "TTree.h"
47 #include "TObjString.h"
48 #include "TSystem.h"
49 
50 #include "reportingUtils.hpp"
51 
52 
53 using namespace std;
54 using namespace rpwa;
55 
56 
57 void
58 usage(const string& progName,
59  const int errCode = 0)
60 {
61  cerr << "calculates difference of two amplitude files an writes result into tree" << endl
62  << endl
63  << "usage:" << endl
64  << progName
65  << " [-o out file -v -h] .amp file 1 .amp file 2" << endl
66  << " where:" << endl
67  << " -o file path to output ROOT file (default: none)" << endl
68  << " -v verbose; print debug output (default: false)" << endl
69  << " -h print help" << endl
70  << endl;
71  exit(errCode);
72 }
73 
74 
75 // tries to extract mass bin boundaries from amplitude file path
76 // assumes standard directory layout:
77 // .../[massBinMin].[massBinMax]/{PSP,}AMPS/*.amp
78 bool
79 massBinFromPath(const string& path,
80  unsigned int& massBinMin,
81  unsigned int& massBinMax)
82 {
83  // strip .amp file
84  size_t pos;;
85  if ((pos = path.find_last_of("/")) == string::npos)
86  return false;
87  string massBinDir = path.substr(0, pos);
88  // strip {PSP,}AMP
89  if ((pos = massBinDir.find_last_of("/")) == string::npos)
90  return false;
91  massBinDir = massBinDir.substr(0, pos);
92  // extract [massBinMin].[massBinMax]
93  if ((pos = massBinDir.find_last_of("/")) == string::npos)
94  return false;
95  massBinDir = massBinDir.substr(pos + 1);
96  // extract mass bins
97  if ((pos = massBinDir.find_last_of(".")) == string::npos)
98  return false;
99  massBinMin = atoi(massBinDir.substr(0, pos ).c_str());
100  massBinMax = atoi(massBinDir.substr(pos + 1).c_str());
101  if ((massBinMin == 0) or (massBinMax == 0))
102  return false;
103  return true;
104 }
105 
106 
107 int
108 main(int argc,
109  char** argv)
110 {
111  printCompilerInfo();
112  printLibraryInfo ();
113  printSvnVersion ();
114  cout << endl;
115 
116  // parse command line options
117  const string progName = argv[0];
118  string outFileName = "";
119  bool debug = false;
120  string ampFileNames[2] = {"", ""};
121  extern char* optarg;
122  extern int optind;
123  int c;
124  while ((c = getopt(argc, argv, "o:vh")) != -1)
125  switch (c) {
126  case 'o':
127  outFileName = optarg;
128  break;
129  case 'v':
130  debug = true;
131  break;
132  case 'h':
133  default:
134  usage(progName);
135  }
136 
137  // get amplitude file names
138  if (argc - optind != 2) {
139  printErr << "you have to specify two amplitude files. aborting." << endl;;
140  usage(progName, 1);
141  }
142  ampFileNames[0] = argv[optind++];
143  ampFileNames[1] = argv[optind++];
144 
145  // open amplitude files
146  printInfo << "comparing amplitude files '" << ampFileNames[0] << "' and "
147  << "'" << ampFileNames[1] << endl;
148  ifstream ampFiles[2];
149  for (unsigned int i = 0; i < 2; ++i) {
150  ampFiles[i].open(ampFileNames[i].c_str());
151  if (not ampFiles[i]) {
152  printErr << "cannot open amplitude file '" << ampFileNames[i] << "'. aborting." << endl;
153  exit(1);
154  }
155  }
156 
157  // read amplitudes into memory
158  vector<complex<double> > amps[2];
159  for (unsigned int i = 0; i < 2; ++i) {
160  complex<double> amp;
161  while (ampFiles[i].read((char*) &amp, sizeof(complex<double>)))
162  amps[i].push_back(amp);
163  }
164  if (amps[0].size() != amps[1].size())
165  printWarn << "the two amplitude files have different number of amplitudes "
166  << "(" << amps[0].size() << " vs. " << amps[1].size() << ")." << endl;
167  const unsigned int nmbAmps = min(amps[0].size(), amps[1].size());
168 
169  // open output file
170  TFile* outFile = 0;
171  if (outFileName != "") {
172  printInfo << "writing difference tree to '" << outFileName << "'" << endl;
173  outFile = TFile::Open(outFileName.c_str(), "RECREATE");
174  }
175 
176  // create tree
177  TTree* tree = new TTree("ampDiffTree", "ampDiffTree");
178  if (not tree) {
179  printErr << "problems creating tree 'ampDiffTree' " << "in file '" << outFileName << "'" << endl;
180  return false;
181  }
182 
183  // create leaf variables
184  const string cmd = "basename '" + ampFileNames[0] + "'";
185  TObjString* ampName = new TObjString(gSystem->GetFromPipe(cmd.c_str()));
186  UInt_t eventNmb;
187  UInt_t massBinMin = 0, massBinMax = 0; // [MeV/c^2]
188  Double_t valReal[2], valImag[2];
189  Double_t absDiffReal, absDiffImag;
190  Double_t relDiffReal, relDiffImag;
191 
192  // get mass bin boundaries from file paths
193  if (not massBinFromPath(ampFileNames[0], massBinMin, massBinMax)
194  and not massBinFromPath(ampFileNames[1], massBinMin, massBinMax))
195  printWarn << "cannot determine mass bin boundaries from file paths" << endl;
196  else
197  printInfo << "extracted mass bin boundaries [" << massBinMin << ", " << massBinMax << "] from file paths" << endl;
198 
199  // connect leaf variables to tree branches
200  const int split = 0;
201  const int bufSize = 256000;
202  tree->Branch("ampName", "TObjString", &ampName, bufSize, split);
203  tree->Branch("eventNmb", &eventNmb, "eventNmb/i");
204  tree->Branch("massBinMin", &massBinMin, "massBinMin/i");
205  tree->Branch("massBinMax", &massBinMax, "massBinMax/i");
206  tree->Branch("valReal", valReal, "valReal[2]/D");
207  tree->Branch("valImag", valImag, "valImag[2]/D");
208  tree->Branch("absDiffReal", &absDiffReal, "absDiffReal/D");
209  tree->Branch("absDiffImag", &absDiffImag, "absDiffImag/D");
210  tree->Branch("relDiffReal", &relDiffReal, "relDiffReal/D");
211  tree->Branch("relDiffImag", &relDiffImag, "relDiffImag/D");
212 
213  // compare amplitudes
214  double maxAbsDiff = 0;
215  long int maxAbsDiffIndex = -1;
216  double maxRelDiff = 0;
217  long int maxRelDiffIndex = -1;
218  for (unsigned int i = 0; i < nmbAmps; ++i) {
219  const complex<double> absDiff = amps[0][i] - amps[1][i];
220  const complex<double> relDiff = complex<double>(absDiff.real() / amps[0][i].real(),
221  absDiff.imag() / amps[0][i].imag());
222  // fill tree
223  if (outFile) {
224  eventNmb = i;
225  valReal[0] = amps[0][i].real();
226  valImag[0] = amps[0][i].imag();
227  valReal[1] = amps[1][i].real();
228  valImag[1] = amps[1][i].imag();
229  absDiffReal = absDiff.real();
230  absDiffImag = absDiff.imag();
231  relDiffReal = relDiff.real();
232  relDiffImag = relDiff.imag();
233  tree->Fill();
234  }
235  // print amplitudes
236  if (debug) {
237  const unsigned int nmbDigits = numeric_limits<double>::digits10 + 1;
238  ostringstream s;
239  s.precision(nmbDigits);
240  s.setf(ios_base::scientific, ios_base::floatfield);
241  s << " " << setw(49) << amps[0][i] << " - " << setw(49) << amps[1][i]
242  << " = " << setw(49) << absDiff << ", relative "
243  << "(" << setw(23) << relDiff.real() << ", " << setw(23) << relDiff.imag() << " )";
244  cout << s.str() << endl;
245  }
246  // compute maximum deviations
247  const double absMax = max(fabs(absDiffReal), fabs(absDiffImag));
248  const double relMax = max(fabs(relDiffReal), fabs(relDiffImag));
249  if (absMax > maxAbsDiff) {
250  maxAbsDiff = absMax;
251  maxAbsDiffIndex = i;
252  }
253  if (relMax > maxRelDiff) {
254  maxRelDiff = relMax;
255  maxRelDiffIndex = i;
256  }
257  }
258  printInfo << "maximum observed deviations: absolute = " << maxPrecision(maxAbsDiff) << " "
259  << "(event " << maxAbsDiffIndex << "), relative = " << maxPrecision(maxRelDiff) << " "
260  << "(event " << maxRelDiffIndex << ")" << endl;
261 
262  if (outFile) {
263  outFile->Write();
264  outFile->Close();
265  delete outFile;
266  printInfo << "wrote difference tree to '" << outFileName << "'" << endl;
267  }
268  return 0;
269 }
Generated by   doxygen 1.8.1.2