ROOTPWA
calcAmplitudes.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 // reads in data files in .evt or ROOT tree format, calculates
29 // amplitudes for each event based on given key file, and writes
30 // out amplitudes in ROOT tree or PWA2000 binary or ascii format
31 //
32 //
33 // Author List:
34 // Boris Grube TUM (original author)
35 //
36 //
37 //-------------------------------------------------------------------------
38 
39 
40 #include <iostream>
41 #include <fstream>
42 #include <unistd.h>
43 #include <vector>
44 #include <complex>
45 
46 #include "TROOT.h"
47 #include "TFile.h"
48 #include "TChain.h"
49 #include "TTree.h"
50 #include "TBranch.h"
51 #include "TClonesArray.h"
52 #include "TStopwatch.h"
53 
54 #include "fileUtils.hpp"
55 #include "particleDataTable.h"
56 #include "evtTreeHelper.h"
57 #include "waveDescription.h"
58 #include "amplitudeTreeLeaf.h"
59 
60 
61 using namespace std;
62 using namespace rpwa;
63 
64 
65 void
66 usage(const string& progName,
67  const int errCode = 0)
68 {
69 
70  cerr << "calculates decay amplitudes for given wave for events in input data files and" << endl
71  << "writes amplitudes to file" << endl
72  << endl
73  << "usage:" << endl
74  << progName
75  << " -k key file [-n max. # of events -p PDG file -o output file -m amplitude leaf name "
76  << "-a -t tree name -l leaf names -r target particle name -v -h] "
77  << "input data file(s) (.evt or .root format)" << endl
78  << " where:" << endl
79  << " -k file path to key file" << endl
80  << " -n # maximum number of events to read (default: all)" << endl
81  << " -p file path to particle data table file (default: ./particleDataTable.txt)" << endl
82  << " -o file path to amplitude file (.amp or .root format; default: ./out.root)" << endl
83  << " -m amplitude leaf name (default: 'amplitude')" << endl
84  << " -a write .amp files in ASCII format (default: binary)" << endl
85  << " -t name name of tree in ROOT data files (default: rootPwaEvtTree)" << endl
86  << " -l names semicolon separated object/leaf names in input data" << endl
87  << " (default: 'prodKinParticles;prodKinMomenta;decayKinParticles;decayKinMomenta')" << endl
88  << " -v verbose; print debug output (default: false)" << endl
89  << " -h print help" << endl
90  << endl;
91  exit(errCode);
92 }
93 
94 
95 int
96 main(int argc,
97  char** argv)
98 {
99  printCompilerInfo();
100  printLibraryInfo ();
101  printSvnVersion ();
102  cout << endl;
103 
104 #ifdef USE_STD_COMPLEX_TREE_LEAFS
105  // force loading predefined std::complex dictionary
106  // see http://root.cern.ch/phpBB3/viewtopic.php?f=5&t=9618&p=50164
107  gROOT->ProcessLine("#include <complex>");
108 #endif
109 
110  // parse command line options
111  const string progName = argv[0];
112  string keyFileName = "";
113  long int maxNmbEvents = -1;
114  string pdgFileName = "./particleDataTable.txt";
115  string ampFileName = "./out.root";
116  string ampLeafName = "amplitude";
117  bool asciiOutput = false;
118  string inTreeName = "rootPwaEvtTree";
119  string leafNames = "prodKinParticles;prodKinMomenta;decayKinParticles;decayKinMomenta";
120  bool newKeyFileNameConvention = false;
121  bool debug = false;
122  const long int treeCacheSize = 1000000; // 1 MByte ROOT tree read cache size
123  extern char* optarg;
124  extern int optind;
125  int c;
126  while ((c = getopt(argc, argv, "k:n:p:o:m:at:l:r:vh")) != -1)
127  switch (c) {
128  case 'k':
129  keyFileName = optarg;
130  break;
131  case 'n':
132  maxNmbEvents = atol(optarg);
133  break;
134  case 'p':
135  pdgFileName = optarg;
136  break;
137  case 'o':
138  ampFileName = optarg;
139  break;
140  case 'm':
141  ampLeafName = optarg;
142  break;
143  case 'a':
144  asciiOutput = true;
145  break;
146  case 't':
147  inTreeName = optarg;
148  break;
149  case 'l':
150  leafNames = optarg;
151  break;
152  case 'v':
153  debug = true;
154  break;
155  case 'h':
156  default:
157  usage(progName);
158  }
159 
160  // get input file names
161  if (optind >= argc) {
162  printErr << "you need to specify at least one data file to process. aborting." << endl;;
163  usage(progName, 1);
164  }
165  vector<string> rootFileNames;
166  vector<string> evtFileNames;
167  while (optind < argc) {
168  const string fileName = argv[optind++];
169  const string fileExt = extensionFromPath(fileName);
170  if (fileExt == "root")
171  rootFileNames.push_back(fileName);
172  else if (fileExt == "evt")
173  evtFileNames.push_back(fileName);
174  else
175  printWarn << "input file '" << fileName << "' is neither a .root nor a .evt file. "
176  << "skipping." << endl;
177  }
178  if ((rootFileNames.size() == 0) and (evtFileNames.size() == 0)) {
179  printErr << "none of the specified input files is a .root or .evt file. aborting.";
180  usage(progName, 1);
181  }
182 
183  // get object and leaf names for event data
184  string prodKinPartNamesObjName, prodKinMomentaLeafName;
185  string decayKinPartNamesObjName, decayKinMomentaLeafName;
186  parseLeafAndObjNames(leafNames, prodKinPartNamesObjName, prodKinMomentaLeafName,
187  decayKinPartNamesObjName, decayKinMomentaLeafName);
188 
189  // open .root and .evt files for event data
190  vector<TTree*> inTrees;
191  TClonesArray* prodKinPartNames = 0;
192  TClonesArray* decayKinPartNames = 0;
193  if (not openRootEvtFiles(inTrees, prodKinPartNames, decayKinPartNames,
194  rootFileNames, evtFileNames,
195  inTreeName, prodKinPartNamesObjName, prodKinMomentaLeafName,
196  decayKinPartNamesObjName, decayKinMomentaLeafName, debug)) {
197  printErr << "problems opening input file(s). aborting." << endl;
198  exit(1);
199  }
200 
201  // initialize particle data table
202  particleDataTable::readFile(pdgFileName);
203 
204  // parse key file and create amplitude instance
205  if (keyFileName == "") {
206  printErr << "no key file specified. aborting." << endl;
207  usage(progName, 1);
208  }
209  waveDescription waveDesc;
210  isobarAmplitudePtr amplitude;
211  if ( not waveDesc.parseKeyFile(keyFileName)
212  or not waveDesc.constructAmplitude(amplitude)) {
213  printErr << "problems constructing decay topology from key file '" << keyFileName << "'. "
214  << "aborting." << endl;
215  exit(1);
216  }
217  printInfo << *amplitude;
218 
219  // create output file for amplitudes
220  bool writeRootFormat = false;
221 #ifdef USE_STD_COMPLEX_TREE_LEAFS
222  const string ampFileExt = extensionFromPath(ampFileName);
223  if (ampFileExt == "root")
224  writeRootFormat = true;
225  else if (ampFileExt == "amp")
226  writeRootFormat = false;
227  else {
228  printErr << "specified amplitude file '" << ampFileName << "' is neither a .root "
229  << "nor a .amp file. aborting." << endl;
230  usage(progName);
231  }
232 #endif
233  printInfo << "creating amplitude file '" << ampFileName << "'";
234  ofstream* ampFilePlain = 0;
235  TFile* ampFileRoot = 0;
236 #ifdef USE_STD_COMPLEX_TREE_LEAFS
237  amplitudeTreeLeaf* ampTreeLeaf = 0;
238  TTree* ampTree = 0;
239  if (writeRootFormat) {
240  cout << endl;
241  ampFileRoot = TFile::Open(ampFileName.c_str(), "RECREATE");
242  // write wave description
243  const string waveName = waveDesc.waveNameFromTopology(*(amplitude->decayTopology()),
244  newKeyFileNameConvention);
245  waveDesc.Write((waveName + ".key").c_str());
246  // create output tree
247  ampTreeLeaf = new amplitudeTreeLeaf();
248  const string ampTreeName = waveName + ".amp";
249  ampTree = new TTree(ampTreeName.c_str(), ampTreeName.c_str());
250  const int splitLevel = 99;
251  const int bufSize = 256000;
252  ampTree->Branch(ampLeafName.c_str(), &ampTreeLeaf, bufSize, splitLevel);
253  } else
254 #endif
255  {
256  cout << "; " << ((asciiOutput) ? "ASCII" : "binary") << " mode" << endl;
257  ampFilePlain = new ofstream(ampFileName.c_str());
258  }
259  if ( (writeRootFormat and not ampFileRoot)
260  or (not writeRootFormat and not ampFilePlain and not *ampFilePlain)) {
261  printErr << "cannot create amplitude file '" << ampFileName << "'. aborting." << endl;
262  exit(1);
263  }
264 
265  // read data from tree(s) and calculate decay amplitudes
266  TStopwatch timer;
267  timer.Reset();
268  timer.Start();
269  vector<complex<double> > ampValues;
270  for (unsigned int i = 0; i < inTrees.size(); ++i) {
271  printInfo << "processing ";
272  if ((rootFileNames.size() > 0) and (i == 0))
273  cout << "chain of .root files";
274  else
275  cout << ".evt tree[" << ((rootFileNames.size() > 0) ? i : i + 1) << "]";
276  cout << endl;
277  if (not processTree(*inTrees[i], *prodKinPartNames, *decayKinPartNames,
278  amplitude, ampValues, maxNmbEvents - ampValues.size(),
279  prodKinMomentaLeafName, decayKinMomentaLeafName))
280  printWarn << "problems reading tree" << endl;
281  }
282  printSucc << "calculated amplitudes for " << ampValues.size() << " events" << endl;
283  timer.Stop();
284  printInfo << "this job consumed: ";
285  timer.Print();
286 
287  // write amplitudes to output file
288  for (unsigned int i = 0; i < ampValues.size(); ++i) {
289 #ifdef USE_STD_COMPLEX_TREE_LEAFS
290  if (ampFileRoot) {
291  ampTreeLeaf->setAmp(ampValues[i]);
292  ampTree->Fill();
293  }
294 #endif
295  if (ampFilePlain) {
296  if (asciiOutput)
297  *ampFilePlain << setprecision(numeric_limits<double>::digits10 + 1) << ampValues[i] << endl;
298  else
299  ampFilePlain->write((char*)(&ampValues[i]), sizeof(complex<double>));
300  }
301  }
302 #ifdef USE_STD_COMPLEX_TREE_LEAFS
303  if (ampFileRoot) {
304  printInfo << "optimizing tree" << endl;
305  //ampTree->Print();
306  ampTree->OptimizeBaskets(treeCacheSize, 1, "d");
307  ampTree->Write();
308  //ampTree->Print();
309  ampFileRoot->Close();
310  delete ampFileRoot;
311  }
312 #endif
313  if (ampFilePlain) {
314  ampFilePlain->close();
315  delete ampFilePlain;
316  }
317  printSucc << "wrote " << ampValues.size() << " amplitude values to "
318  << "'" << ampFileName << "'";
319  if (not writeRootFormat)
320  cout << " in " << ((asciiOutput) ? "ASCII" : "binary")
321  << " mode";
322  cout << endl;
323 
324  // clean up
325  for (unsigned int i = 0; i < inTrees.size(); ++i)
326  delete inTrees[i];
327  inTrees.clear();
328 
329  return 0;
330 }
Generated by   doxygen 1.8.1.2