ROOTPWA
evtTreeHelper.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 // helper functions that convert between standard ASCII PWA2000
29 // .evt files and the new ROOT tree format
30 //
31 //
32 // Author List:
33 // Boris Grube TUM (original author)
34 //
35 //
36 //-------------------------------------------------------------------------
37 
38 
39 #include <fstream>
40 #include <sstream>
41 #include <string>
42 #include <cassert>
43 #include <algorithm>
44 #include <map>
45 
46 #include <boost/tokenizer.hpp>
47 #include <boost/progress.hpp>
48 #include <boost/bimap.hpp>
49 #include <boost/assign/list_inserter.hpp>
50 
51 #include "TFile.h"
52 #include "TTree.h"
53 #include "TTreePerfStats.h"
54 #include "TChain.h"
55 #include "TClonesArray.h"
56 #include "TObjString.h"
57 #include "TVector3.h"
58 
59 #include "reportingUtilsRoot.hpp"
60 #include "conversionUtils.hpp"
61 #include "particleDataTable.h"
62 #include "isobarDecayTopology.h"
63 #include "isobarHelicityAmplitude.h"
64 #include "evtTreeHelper.h"
65 
66 
67 using namespace std;
68 using namespace boost;
69 using namespace boost::bimaps;
70 
71 
72 namespace rpwa {
73 
74 
75  bool
76  checkParticleCharge(const size_t lineNmb,
77  const int id,
78  const string& name,
79  const int chargeToCheck)
80  {
81  if (name == "unknown") {
82  printWarn << "error reading data line " << lineNmb << ": unknown particle" << endl;
83  return false;
84  } else {
85  // check charge
86  int charge;
87  particleProperties::chargeFromName(name, charge);
88  if (chargeToCheck != charge) {
89  printWarn << "error reading data line " << lineNmb << ": "
90  << "GEANT particle ID of " << id << " corresponds to particle "
91  << "'" << name << "'. this inconsistent with the charge of "
92  << "'" << chargeToCheck << "' in data file.";
93  return false;
94  }
95  }
96  return true;
97  }
98 
99 
100  double
101  getParticleMass(const string& name)
102  {
103  rpwa::particleDataTable& pdt = rpwa::particleDataTable::instance();
104  const particleProperties* prop = 0;
105  if (pdt.isInTable(name))
106  prop = pdt.entry(name);
107  else {
108  const string n = particleProperties::stripChargeFromName(name);
109  if (pdt.isInTable(n))
110  prop = pdt.entry(n);
111  }
112  if (not prop) {
113  printWarn << "neither particle '" << name << "' "
114  << "nor '" << particleProperties::stripChargeFromName(name) << "' "
115  << "are in particle data table. using mass 0." << endl;
116  return 0;
117  }
118  return prop->mass();
119  }
120 
121 
122  void
123  parseLeafAndObjNames(const string& cmdLineString,
124  string& prodKinPartNamesObjName,
125  string& prodKinMomentaLeafName,
126  string& decayKinPartNamesObjName,
127  string& decayKinMomentaLeafName)
128  {
129  typedef tokenizer<char_separator<char> > tokenizer;
130  char_separator<char> separator(";");
131  tokenizer nameTokens(cmdLineString, separator);
132  tokenizer::iterator nameToken = nameTokens.begin();
133  prodKinPartNamesObjName = *nameToken;
134  prodKinMomentaLeafName = *(++nameToken);
135  decayKinPartNamesObjName = *(++nameToken);
136  decayKinMomentaLeafName = *(++nameToken);
137  printInfo << "using the following object/leaf names:" << endl
138  << " production kinematics: "
139  << "particle names = '" << prodKinPartNamesObjName << "', "
140  << "momenta = '" << prodKinMomentaLeafName << "'" << endl
141  << " decay kinematics: "
142  << "particle names = '" << decayKinPartNamesObjName << "', "
143  << "momenta = '" << decayKinMomentaLeafName << "'" << endl;
144  }
145 
146 
147  bool
148  getParticleNamesFromRootFile(TFile& inFile,
149  TClonesArray*& prodKinPartNames, // array of particle names to be filled
150  TClonesArray*& decayKinPartNames, // array of particle names to be filled
151  const string& inTreeName,
152  const string& prodKinPartNamesObjName,
153  const string& decayKinPartNamesObjName)
154  {
155  // read particle names from first root file
156  inFile.GetObject(prodKinPartNamesObjName.c_str(), prodKinPartNames );
157  inFile.GetObject(decayKinPartNamesObjName.c_str(), decayKinPartNames);
158  if (prodKinPartNames and decayKinPartNames) {
159  TVector3::Class()->IgnoreTObjectStreamer(true);
160  return true;
161  }
162  // if particle names are not in file try first tree entry for backward compatibility
163  TTree* inTree = 0;
164  inFile.GetObject(inTreeName.c_str(), inTree);
165  if (not prodKinPartNames) {
166  TBranch* prodKinPartNamesBr = 0;
167  if (inTree->SetBranchAddress(prodKinPartNamesObjName.c_str(),
168  &prodKinPartNames, &prodKinPartNamesBr) >= 0)
169  if (inTree->GetEntries() > 0)
170  if (inTree->LoadTree(0) >= 0)
171  prodKinPartNamesBr->GetEntry(0);
172  }
173  if (not decayKinPartNames) {
174  TBranch* decayKinPartNamesBr = 0;
175  if (inTree->SetBranchAddress(decayKinPartNamesObjName.c_str(),
176  &decayKinPartNames, &decayKinPartNamesBr) >= 0)
177  if (inTree->GetEntries() > 0)
178  if (inTree->LoadTree(0) >= 0)
179  decayKinPartNamesBr->GetEntry(0);
180  }
181  if (prodKinPartNames and decayKinPartNames)
182  return true;
183  else
184  return false;
185  }
186 
187 
188  bool
189  openRootEvtFiles(vector<TTree*>& inTrees, // array of trees from .root and .evt files
190  TClonesArray*& prodKinPartNames, // array of particle names to be filled
191  TClonesArray*& decayKinPartNames, // array of particle names to be filled
192  const vector<string>& rootFileNames, // .root files to be opened
193  const vector<string>& evtFileNames, // .evt files to be converted to trees
194  const string& inTreeName,
195  const string& prodKinPartNamesObjName,
196  const string& prodKinMomentaLeafName,
197  const string& decayKinPartNamesObjName,
198  const string& decayKinMomentaLeafName,
199  const bool debug)
200  {
201  // open root files and build chain
202  TChain* inChain = 0;
203  if (rootFileNames.size() > 0) {
204  inChain = new TChain(inTreeName.c_str());
205  for (unsigned int i = 0; i < rootFileNames.size(); ++i) {
206  printInfo << "opening ROOT input file '" << rootFileNames[i] << "'" << endl;
207  if (inChain->Add(rootFileNames[i].c_str()) < 1)
208  printWarn << "no events in ROOT input file '" << rootFileNames[i] << "'" << endl;
209  }
210  // read particle names from first root file
211  TFile* inFile = inChain->GetFile(); // opens first file
212  if (not inFile)
213  printWarn << "could not open file '" << rootFileNames[0] << "'" << endl;
214  else {
215  if (not getParticleNamesFromRootFile(*inFile, prodKinPartNames, decayKinPartNames,
216  inTreeName,
217  prodKinPartNamesObjName, decayKinPartNamesObjName)) {
218  printErr << "cannot find production and/or decay kinematics particle names "
219  << "in input file '" << inChain->GetFile()->GetName() << "'" << endl;
220  return false;
221  }
222  }
223  }
224  if (inChain)
225  inTrees.push_back(inChain);
226 
227  // convert .evt files to root trees
228  for (unsigned int i = 0; i < evtFileNames.size(); ++i) {
229  printInfo << "opening .evt input file '" << evtFileNames[i] << "'" << endl;
230  ifstream evtFile(evtFileNames[i].c_str());
231  if (not evtFile or not evtFile.good()) {
232  printWarn << "cannot open .evt input file '" << evtFileNames[i] << "'. skipping." << endl;
233  continue;
234  }
235  printInfo << "converting .evt input file '" << evtFileNames[i] << "' "
236  << "into memory resident tree. this might reduce performance. "
237  << "ROOT input format is recommended." << endl;
238  // create tree
239  TTree* tree = new TTree(inTreeName.c_str(), inTreeName.c_str());
240  if (not tree) {
241  printErr << "problems creating tree '" << inTreeName << "'. skipping." << endl;
242  continue;
243  }
244  TClonesArray prodNames ("TObjString"); // names of production kinematics particles read from .evt file
245  TClonesArray decayNames("TObjString"); // names of decay kinematics particles read from .evt file
246  if (fillTreeFromEvt(evtFile, *tree, prodNames, decayNames, -1,
247  prodKinMomentaLeafName, decayKinMomentaLeafName, debug))
248  inTrees.push_back(tree);
249  else {
250  printWarn << "problems creating tree from .evt input file '" << evtFileNames[i] << "' "
251  << "skipping." << endl;
252  }
253 
254  // set particle names
255  if (prodKinPartNames) {
256  // check consistency
257  int j;
258  if (prodNames.GetEntriesFast() != prodKinPartNames->GetEntriesFast())
259  goto prodKinPartNamesInconsistent;
260  for (j = 0; j < prodNames.GetEntriesFast(); ++j) {
261  if (not (((TObjString*)prodNames[j])->GetString()
262  == ((TObjString*)(*prodKinPartNames)[j])->GetString()))
263  goto prodKinPartNamesInconsistent;
264  prodKinPartNamesInconsistent:
265  {
266  printErr << "production kinematics particle names read from input files "
267  << "are inconsistent:" << endl
268  << " expected:";
269  for (int k = 0; k < prodKinPartNames->GetEntriesFast(); ++k)
270  cout << " " << ((TObjString*)(*prodKinPartNames)[k])->GetString().Data()
271  << "[" << k << "]";
272  cout << endl << " read from '" << evtFileNames[i] << "':";
273  for (int k = 0; k < prodNames.GetEntriesFast(); ++k)
274  cout << " " << ((TObjString*)prodNames[k])->GetString().Data() << "[" << k << "]";
275  cout << endl << " check input files." << endl;
276  return false;
277  }
278  }
279  } else {
280  prodKinPartNames = new TClonesArray("TObjString");
281  *prodKinPartNames = prodNames;
282  }
283  if (decayKinPartNames) {
284  // check consistency
285  int j;
286  if (decayNames.GetEntriesFast() != decayKinPartNames->GetEntriesFast())
287  goto decayKinPartNamesInconsistent;
288  for (j = 0; j < decayNames.GetEntriesFast(); ++j) {
289  if (not (((TObjString*)decayNames[j])->GetString()
290  == (((TObjString*)(*decayKinPartNames)[j])->GetString())))
291  goto decayKinPartNamesInconsistent;
292  decayKinPartNamesInconsistent:
293  {
294  printErr << "decay kinematics particle names read from input files "
295  << "are inconsistent:" << endl
296  << " expected:";
297  for (int k = 0; k < decayKinPartNames->GetEntriesFast(); ++k)
298  cout << " " << ((TObjString*)(*decayKinPartNames)[k])->GetString().Data()
299  << "[" << k << "]";
300  cout << endl << " read from '" << evtFileNames[i] << "':";
301  for (int k = 0; k < decayNames.GetEntriesFast(); ++k)
302  cout << " " << ((TObjString*)decayNames[k])->GetString().Data() << "[" << k << "]";
303  cout << endl << " check input files." << endl;
304  return false;
305  }
306  }
307  } else {
308  decayKinPartNames = new TClonesArray("TObjString");
309  *decayKinPartNames = decayNames;
310  }
311  } // loop over .evt files
312 
313  // make sure particle names are specified
314  bool success = true;
315  if (not prodKinPartNames) {
316  printErr << "no production kinematics particle names were found in input files." << endl;
317  success = false;
318  }
319  if (not decayKinPartNames) {
320  printWarn << "no decay kinematics particle names were found in input files." << endl;
321  success = false;
322  }
323 
324  return success;
325  }
326 
327 
328  bool
329  fillTreeFromEvt(istream& inEvt,
330  TTree& outTree, // tree to be filled
331  TClonesArray& prodKinPartNames, // array of particle names to be filled
332  TClonesArray& decayKinPartNames, // array of particle names to be filled
333  const long int maxNmbEvents,
334  const string& prodKinMomentaLeafName,
335  const string& decayKinMomentaLeafName,
336  const bool debug,
337  const long int treeCacheSize)
338  {
339  if (not inEvt or not inEvt.good()) {
340  printWarn << "cannot read from input stream" << endl;
341  return false;
342  }
343 
344  // create leaf variables
345  TVector3::Class()->IgnoreTObjectStreamer(true);
346  TClonesArray* prodKinMomenta = new TClonesArray("TVector3");
347  TClonesArray* decayKinMomenta = new TClonesArray("TVector3");
348 
349  // connect leaf variables to tree branches
350  const int splitLevel = 99;
351  const int bufSize = 256000;
352  outTree.Branch(prodKinMomentaLeafName.c_str(), "TClonesArray", &prodKinMomenta, bufSize, splitLevel);
353  outTree.Branch(decayKinMomentaLeafName.c_str(), "TClonesArray", &decayKinMomenta, bufSize, splitLevel);
354 
355  // loop over events and fill tree
356  vector<string> prodNames; // names of production kinematics particles read from .evt file
357  vector<string> decayNames; // names of decay kinematics particles read from .evt file
358  bool success = true;
359  long int countEvents = 0;
360  long int countLines = 0;
361  inEvt.seekg(0, ios::end);
362  long int fileLength = inEvt.tellg();
363  inEvt.seekg(0, ios::beg);
364  progress_display* progressIndicator = (not debug) ? new progress_display(fileLength, cout, "") : 0;
365  streampos lastPos = inEvt.tellg();
366  while (inEvt.good()) {
367  string line;
368 
369  // read number of particles
370  int nmbParticles = 0;
371  if (getline(inEvt, line)) {
372  ++countLines;
373  stringstream lineStream(line);
374  int n;
375  if (lineStream >> n)
376  nmbParticles = n;
377  else {
378  printWarn << "event " << countEvents + 1 << ": error reading number of particles "
379  << "from line " << countLines << ": " << line << endl;
380  success = false;
381  }
382  } else
383  break;
384  assert(nmbParticles > 0);
385  if (debug)
386  printDebug << "# of particles = " << nmbParticles << endl;
387 
388  // read production kinematics data (beam + fixed target)
389  prodNames.clear();
390  prodKinMomenta->Clear();
391  if (getline(inEvt, line)) {
392  ++countLines;
393  stringstream lineStream(line);
394  int id = 0, charge = 0;
395  double momX = 0, momY = 0, momZ = 0, E = 0;
396  if (lineStream >> id >> charge >> momX >> momY >> momZ >> E) {
397  const string partName = particleDataTable::particleNameFromGeantId(id);
398  prodNames.push_back(partName);
399  if (not checkParticleCharge(countLines, id, partName, charge))
400  success = false;
401  new((*prodKinMomenta)[0]) TVector3(momX, momY, momZ);
402  } else {
403  printWarn << "event " << countEvents + 1 << ": error reading beam data "
404  << "from line " << countLines << ": " << line << endl;
405  success = false;
406  }
407  } else
408  break;
409 
410  // check consistency
411  const int nmbProdKinPart = prodNames.size();
412  assert((nmbProdKinPart > 0) and (nmbProdKinPart == prodKinMomenta->GetEntriesFast()));
413  if (debug) {
414  printDebug << nmbProdKinPart << " production kinematics particles:" << endl;
415  for (int i = 0; i < nmbProdKinPart; ++i)
416  cout << " particle[" << i << "]: "
417  << prodNames[i] << "; " << *((TVector3*)(*prodKinMomenta)[i]) << endl;
418  }
419 
420  // read decay kinematics data
421  decayNames.clear();
422  decayNames.resize(nmbParticles - 1, "");
423  decayKinMomenta->Clear();
424  for (int i = 0; i < nmbParticles - 1; ++i) {
425  if (getline(inEvt, line)) {
426  ++countLines;
427  stringstream lineStream(line);
428  int id, charge;
429  double momX, momY, momZ, E;
430  if (lineStream >> id >> charge >> momX >> momY >> momZ >> E) {
431  const string partName = particleDataTable::particleNameFromGeantId(id);
432  decayNames[i] = partName;
433  if (not checkParticleCharge(countLines, id, partName, charge))
434  success = false;
435  new((*decayKinMomenta)[i]) TVector3(momX, momY, momZ);
436  } else {
437  printWarn << "event " << countEvents + 1 << ": error reading decay kinematics "
438  << "particle[" << i << "] data from line " << countLines << ": " << line << endl;
439  success = false;
440  }
441  } else
442  break;
443  }
444 
445  // check consistency
446  const int nmbDecayKinPart = decayNames.size();
447  assert((nmbDecayKinPart > 0) and (nmbDecayKinPart == decayKinMomenta->GetEntriesFast()));
448  if (debug) {
449  printDebug << nmbDecayKinPart << " decay kinematics particles:" << endl;
450  for (int i = 0; i < nmbDecayKinPart; ++i)
451  cout << " particle[" << i << "]: "
452  << decayNames[i] << "; " << *((TVector3*) (*decayKinMomenta)[i]) << endl;
453  }
454 
455  outTree.Fill();
456  if (countEvents == 0) {
457  for (unsigned int i = 0; i < prodNames.size(); ++i)
458  new(prodKinPartNames[i]) TObjString(prodNames[i].c_str());
459  for (unsigned int i = 0; i < decayNames.size(); ++i)
460  new(decayKinPartNames[i]) TObjString(decayNames[i].c_str());
461  } else {
462  // check that particle lists do not change
463  const int nmbIsPart = prodNames.size();
464  const int nmbIsPartExpected = prodKinPartNames.GetEntriesFast();
465  const int nmbFsPart = decayNames.size();
466  const int nmbFsPartExpected = decayKinPartNames.GetEntriesFast ();
467  if ((nmbIsPart != nmbIsPartExpected) || (nmbFsPart != nmbFsPartExpected)) {
468  printErr << "number of particles changed in event " << countEvents << ". "
469  << "producion kinematics: " << nmbIsPart << " particles "
470  << "(" << nmbIsPartExpected << " expected); "
471  << "decay kinematics: " << nmbFsPart << " particles "
472  << "(" << nmbFsPartExpected << " expected). "
473  << "check .evt file. aborting." << endl;
474  throw;
475  }
476  for (int i = 0; i < nmbIsPart; ++i) {
477  const string name = prodNames[i];
478  const string nameExpected = ((TObjString*)prodKinPartNames[i])->GetString().Data();
479  if (name != nameExpected) {
480  printErr << "particle[" << i << "] name mismatch in production kinematics in event "
481  << countEvents << ": '" << name << "', expected '" << nameExpected << "'. "
482  << "check .evt file. aborting." << endl;
483  throw;
484  }
485  }
486  for (int i = 0; i < nmbFsPart; ++i) {
487  const string name = decayNames[i];
488  const string nameExpected = ((TObjString*)decayKinPartNames[i])->GetString().Data();
489  if (name != nameExpected) {
490  printErr << "particle[" << i << "] name mismatch in decay kinematics in event "
491  << countEvents << ": '" << name << "', expected '" << nameExpected << "'. "
492  << "check .evt file. aborting." << endl;
493  throw;
494  }
495  }
496  }
497 
498  ++countEvents;
499  if (progressIndicator)
500  (*progressIndicator) += inEvt.tellg() - lastPos;
501  lastPos = inEvt.tellg();
502  if ((maxNmbEvents > 0) and (countEvents >= maxNmbEvents))
503  break;
504  }
505 
506  printInfo << "optimizing tree" << endl;
507  //outTree.Print();
508  outTree.OptimizeBaskets(treeCacheSize, 1, "d");
509  //outTree.Print();
510 
511  printInfo << "read " << countLines << " lines from input stream and wrote "
512  << countEvents << " events to tree '" << outTree.GetName() << "' "
513  << "assuming fixed target" << endl;
514  return success;
515  }
516 
517 
518  bool
519  writeEvtFromTree(TTree& inTree,
520  ostream& outEvt,
521  const TClonesArray& prodKinPartNames,
522  const TClonesArray& decayKinPartNames,
523  const long int maxNmbEvents,
524  const string& inTreeName,
525  const string& prodKinMomentaLeafName,
526  const string& decayKinMomentaLeafName,
527  const bool debug)
528  {
529  const long int nmbEventsTree = inTree.GetEntries();
530  if (not outEvt) {
531  printWarn << "cannot write to output stream" << endl;
532  return false;
533  }
534  string partClassName = prodKinPartNames.GetClass()->GetName();
535  if (partClassName != "TObjString") {
536  printWarn << "production kinematics particle names are of type '" << partClassName
537  << "' and not TObjString." << endl;
538  return false;
539  }
540  partClassName = decayKinPartNames.GetClass()->GetName();
541  if (partClassName != "TObjString") {
542  printWarn << "decay kinematics particle names are of type '" << partClassName
543  << "' and not TObjString." << endl;
544  return false;
545  }
546 
547  // create leaf variables
548  TClonesArray* prodKinMomenta = 0;
549  TClonesArray* decayKinMomenta = 0;
550 
551  // connect leaf variables to tree branches
552  inTree.SetBranchAddress(prodKinMomentaLeafName.c_str(), &prodKinMomenta );
553  inTree.SetBranchAddress(decayKinMomentaLeafName.c_str(), &decayKinMomenta);
554 
555  // loop over events
556  const long int nmbEvents = ((maxNmbEvents > 0) ? min(maxNmbEvents, nmbEventsTree)
557  : nmbEventsTree);
558  progress_display* progressIndicator = (not debug) ? new progress_display(nmbEvents, cout, "") : 0;
559  for (long int eventIndex = 0; eventIndex < nmbEvents; ++eventIndex) {
560  if (progressIndicator)
561  ++(*progressIndicator);
562 
563  if (inTree.LoadTree(eventIndex) < 0)
564  break;
565  inTree.GetEntry(eventIndex);
566 
567  assert(prodKinMomenta );
568  assert(decayKinMomenta);
569  const int nmbProdKinPart = prodKinPartNames.GetEntriesFast ();
570  const int nmbDecayKinPart = decayKinPartNames.GetEntriesFast();
571  assert(nmbProdKinPart == prodKinMomenta->GetEntriesFast ());
572  assert(nmbDecayKinPart == decayKinMomenta->GetEntriesFast());
573  if (nmbProdKinPart < 1) {
574  printWarn << "array of production kinematics particles does not have any entries. "
575  << "at least entry for beam (index 0) is required. skipping event." << endl;
576  continue;
577  }
578  if (nmbDecayKinPart < 1) {
579  printWarn << "array of decay kinematics particles does not have any entries. "
580  << "at least one final state particle is required. skipping event." << endl;
581  continue;
582  }
583 
584  // write total number of particles
585  outEvt << 1 + nmbDecayKinPart << endl; // PWA2000 supports only beam in production kinematics
586 
587  // write production kinematics (beam only)
588  if (debug)
589  printDebug << "event[" << eventIndex << "]: " << nmbProdKinPart
590  << " production kinematics particles:" << endl;
591  { // only beam
592  assert(prodKinPartNames [0]);
593  assert((*prodKinMomenta)[0]);
594  const string name = ((TObjString*)prodKinPartNames[0])->GetString().Data();
595  const TVector3* mom = dynamic_cast<TVector3*>((*prodKinMomenta)[0]);
596  assert(mom);
597  const double mass = getParticleMass(name);
598  int id, charge;
599  particleDataTable::geantIdAndChargeFromParticleName(name, id, charge);
600  outEvt << setprecision(numeric_limits<double>::digits10 + 1)
601  << id << " " << charge << " " << mom->X() << " " << mom->Y() << " " << mom->Z() << " "
602  << sqrt(mass * mass + mom->Mag2()) << endl;
603  if (debug) {
604  cout << " particle[" << 0 << "]: " << name << ", id = " << id << ", "
605  << "charge = " << charge << "; " << *mom << endl;
606  }
607  }
608 
609  // write decay kinematics
610  if (debug)
611  printDebug << "event[" << eventIndex << "]: " << nmbDecayKinPart
612  << " decay kinematics particles:" << endl;
613  for (unsigned int i = 0; i < (unsigned int)nmbDecayKinPart; ++i) {
614  assert(decayKinPartNames [i]);
615  assert((*decayKinMomenta)[i]);
616  const string name = ((TObjString*)decayKinPartNames[i])->GetString().Data();
617  const TVector3* mom = dynamic_cast<TVector3*>((*decayKinMomenta)[i]);
618  assert(mom);
619  const double mass = getParticleMass(name);
620  int id, charge;
621  particleDataTable::geantIdAndChargeFromParticleName(name, id, charge);
622  outEvt << setprecision(numeric_limits<double>::digits10 + 1)
623  << id << " " << charge << " " << mom->X() << " " << mom->Y() << " " << mom->Z() << " "
624  << sqrt(mass * mass + mom->Mag2()) << endl;
625  if (debug) {
626  cout << " particle[" << i << "]: " << name << ", id = " << id << ", "
627  << "charge = " << charge << "; " << *mom << endl;
628  }
629  }
630 
631  if (debug)
632  cout << endl;
633  }
634 
635  printInfo << "wrote " << nmbEvents << " events to output stream" << endl;
636  return true;
637  }
638 
639 
640  bool
641  processTree(TTree& tree,
642  const TClonesArray& prodKinPartNames,
643  const TClonesArray& decayKinPartNames,
644  const isobarAmplitudePtr& amplitude,
645  vector<complex<double> >& ampValues,
646  const long int maxNmbEvents,
647  const string& prodKinMomentaLeafName,
648  const string& decayKinMomentaLeafName,
649  const bool printProgress,
650  const string& treePerfStatOutFileName,
651  const long int treeCacheSize)
652  {
653  if (not amplitude) {
654  printWarn << "null pointer to isobar decay amplitude. cannot process tree." << endl;
655  return false;
656  }
657  // initialize amplitude
658  amplitude->init();
659  const isobarDecayTopologyPtr& decayTopo = amplitude->decayTopology();
660 
661  // create branch pointers and leaf variables
662  TBranch* prodKinMomentaBr = 0;
663  TBranch* decayKinMomentaBr = 0;
664  TClonesArray* prodKinMomenta = 0;
665  TClonesArray* decayKinMomenta = 0;
666 
667  // connect leaf variables to tree branches
668  tree.SetBranchAddress(prodKinMomentaLeafName.c_str(), &prodKinMomenta, &prodKinMomentaBr );
669  tree.SetBranchAddress(decayKinMomentaLeafName.c_str(), &decayKinMomenta, &decayKinMomentaBr);
670  tree.SetCacheSize(treeCacheSize);
671  tree.AddBranchToCache(prodKinMomentaLeafName.c_str(), true);
672  tree.AddBranchToCache(decayKinMomentaLeafName.c_str(), true);
673  tree.StopCacheLearningPhase();
674  TTreePerfStats* treePerfStats = 0;
675  if (treePerfStatOutFileName != "")
676  treePerfStats = new TTreePerfStats("ioPerf", &tree);
677 
678  // loop over events
679  if (not decayTopo->initKinematicsData(prodKinPartNames, decayKinPartNames)) {
680  printWarn << "problems initializing input data. cannot read input data." << endl;
681  return false;
682  }
683  const long int nmbEventsTree = tree.GetEntries();
684  const long int nmbEvents = ((maxNmbEvents > 0) ? min(maxNmbEvents, nmbEventsTree)
685  : nmbEventsTree);
686  bool success = true;
687  progress_display* progressIndicator = (printProgress) ? new progress_display(nmbEvents, cout, "") : 0;
688  for (long int eventIndex = 0; eventIndex < nmbEvents; ++eventIndex) {
689  if (progressIndicator)
690  ++(*progressIndicator);
691 
692  if (tree.LoadTree(eventIndex) < 0)
693  break;
694  // read only required branches
695  prodKinMomentaBr->GetEntry (eventIndex);
696  decayKinMomentaBr->GetEntry(eventIndex);
697 
698  if (not prodKinMomenta or not decayKinMomenta) {
699  printWarn << "at least one of the input data arrays is a null pointer: "
700  << " production kinematics: " << "momenta = " << prodKinMomenta << endl
701  << " decay kinematics: " << "momenta = " << decayKinMomenta << endl
702  << "skipping event." << endl;
703  success = false;
704  continue;
705  }
706 
707  if (decayTopo->readKinematicsData(*prodKinMomenta, *decayKinMomenta))
708  ampValues.push_back((*amplitude)());
709  else {
710  printWarn << "problems reading event[" << eventIndex << "]" << endl;
711  success = false;
712  }
713  }
714 
715  if (printProgress)
716  tree.PrintCacheStats();
717  if (treePerfStats) {
718  treePerfStats->SaveAs(treePerfStatOutFileName.c_str());
719  delete treePerfStats;
720  }
721  return success;
722  }
723 
724 
725 } // namespace rpwa
Generated by   doxygen 1.8.1.2