ROOTPWA
createDiffTree.C
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 // creates tree with absolute or relative differences of two trees
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <string>
39 #include <cassert>
40 #include <algorithm>
41 
42 #include <boost/progress.hpp>
43 
44 #include "TChain.h"
45 #include "TFile.h"
46 #include "TTree.h"
47 #include "TClonesArray.h"
48 #include "TObjString.h"
49 #include "TVector3.h"
50 
51 #include "reportingUtils.hpp"
52 #include "reportingUtilsRoot.hpp"
53 #include "evtTreeHelper.h"
54 
55 
56 using namespace std;
57 using namespace boost;
58 using namespace rpwa;
59 
60 
61 void
62 fillDiffLeafs(const TClonesArray& partNames,
63  TClonesArray* partMomenta[3],
64  const unsigned int index,
65  const bool absoluteDiff,
66  const bool debug = false)
67 {
68  const TVector3* momenta[2] = {((TVector3*)(*(partMomenta[0]))[index]),
69  ((TVector3*)(*(partMomenta[1]))[index])};
70  TVector3 diff;
71  if (absoluteDiff)
72  diff = *(momenta[0]) - *(momenta[1]);
73  else {
74  diff.SetX((momenta[0]->X() - momenta[1]->X()) / momenta[0]->X());
75  diff.SetY((momenta[0]->Y() - momenta[1]->Y()) / momenta[0]->Y());
76  diff.SetZ((momenta[0]->Z() - momenta[1]->Z()) / momenta[0]->Z());
77  }
78  new((*(partMomenta[2]))[index]) TVector3(diff);
79 
80  if (debug) {
81  const TObjString name = *((TObjString*)partNames[index]);
82  printDebug << "particle[" << index << "]: '" << name.GetString().Data() << "', "
83  << *(momenta[0]) << " vs. " << *(momenta[1]) << ", "
84  << ((absoluteDiff) ? "absolute" : "relative" ) << " difference " << diff << endl;
85  }
86 }
87 
88 
89 bool
90 createDiffTree(const string& inFileNamePatternA = "testEvents.root",
91  const string& inFileNamePatternB = "testEvents2.root",
92  const string& outFileName = "testDiffTree.root",
93  const bool absoluteDiff = true,
94  const long int maxNmbEvents = -1,
95  const string& inTreeName = "rootPwaEvtTree",
96  const string& prodKinPartNamesObjName = "prodKinParticles",
97  const string& prodKinMomentaLeafName = "prodKinMomenta",
98  const string& decayKinPartNamesObjName = "decayKinParticles",
99  const string& decayKinMomentaLeafName = "decayKinMomenta",
100  const bool debug = false,
101  const long int treeCacheSize = 25000000) // 25 MByte ROOT tree read cache
102 {
103  // open input files
104  const string inFileNamePatterns[2] = {inFileNamePatternA, inFileNamePatternB};
105  TTree* inTrees [2] = {0, 0};
106  TClonesArray* prodKinPartNames [2] = {0, 0};
107  TClonesArray* decayKinPartNames [2] = {0, 0};
108  for (unsigned int i = 0; i < 2; ++i) {
109  vector<string> rootFileNames;
110  rootFileNames.push_back(inFileNamePatterns[i]);
111  vector<string> evtFileNames;
112  vector<TTree*> trees;
113  if (not openRootEvtFiles(trees, prodKinPartNames[i], decayKinPartNames[i],
114  rootFileNames, evtFileNames,
115  inTreeName, prodKinPartNamesObjName, prodKinMomentaLeafName,
116  decayKinPartNamesObjName, decayKinMomentaLeafName, debug)) {
117  printErr << "problems opening input file(s). exiting." << endl;
118  return false;
119  }
120  inTrees[i] = trees[0];
121  }
122 
123  // check compatibility of input data
124  const long int nmbEventsInTrees[2] = {inTrees[0]->GetEntries(), inTrees[1]->GetEntries()};
125  if (nmbEventsInTrees[0] != nmbEventsInTrees[1])
126  printWarn << "trees have different number of events: " << nmbEventsInTrees[0] << " vs. "
127  << nmbEventsInTrees[1] << endl;
128  const int nmbProdKinPart [2] = {prodKinPartNames [0]->GetEntriesFast(),
129  prodKinPartNames [1]->GetEntriesFast()};
130  const int nmbDecayKinPart[2] = {decayKinPartNames[0]->GetEntriesFast(),
131  decayKinPartNames[1]->GetEntriesFast()};
132  assert(nmbProdKinPart [0] == nmbProdKinPart [1]);
133  assert(nmbDecayKinPart[0] == nmbDecayKinPart[1]);
134  for (int i = 0; i < nmbProdKinPart[0]; i++)
135  assert( ((TObjString*)(*(prodKinPartNames[0]))[i])->GetString()
136  == ((TObjString*)(*(prodKinPartNames[1]))[i])->GetString());
137  for (int i = 0; i < nmbDecayKinPart[0]; i++)
138  assert( ((TObjString*)(*(decayKinPartNames[0]))[i])->GetString()
139  == ((TObjString*)(*(decayKinPartNames[1]))[i])->GetString());
140 
141  // create output file
142  printInfo << "creating output file '" << outFileName << "'" << endl;
143  TFile* outFile = TFile::Open(outFileName.c_str(), "RECREATE");
144  if (not outFile) {
145  printWarn << "cannot open output file '" << outFileName << "'" << endl;
146  return false;
147  }
148 
149  // copy particle names to output file
150  prodKinPartNames [0]->Write(prodKinPartNamesObjName.c_str (), TObject::kSingleKey);
151  decayKinPartNames[0]->Write(decayKinPartNamesObjName.c_str(), TObject::kSingleKey);
152 
153  // create output tree
154  TTree* tree = new TTree(inTreeName.c_str(), inTreeName.c_str());
155  if (not tree) {
156  printErr << "problems creating tree '" << inTreeName << "' "
157  << "in file '" << outFileName << "'" << endl;
158  return false;
159  }
160 
161  // create leaf variables
162  TClonesArray* prodKinMomenta [3] = {0, 0, new TClonesArray("TVector3", 0)};
163  TClonesArray* decayKinMomenta[3] = {0, 0, new TClonesArray("TVector3", 0)};
164 
165  // connect leaf variables to input tree branches
166  for (unsigned int i = 0; i < 2; ++i) {
167  inTrees[i]->SetBranchAddress(prodKinMomentaLeafName.c_str(), &prodKinMomenta [i]);
168  inTrees[i]->SetBranchAddress(decayKinMomentaLeafName.c_str(), &decayKinMomenta[i]);
169  }
170 
171  // connect leaf variables to output tree branches
172  const int splitLevel = 99;
173  const int bufSize = 256000;
174  tree->Branch(prodKinMomentaLeafName.c_str(), "TClonesArray", &prodKinMomenta [2], bufSize, splitLevel);
175  tree->Branch(decayKinMomentaLeafName.c_str(), "TClonesArray", &decayKinMomenta[2], bufSize, splitLevel);
176 
177  // loop over events
178  long int nmbEvents = min(nmbEventsInTrees[0], nmbEventsInTrees[1]);
179  if (maxNmbEvents > 0)
180  nmbEvents = min(maxNmbEvents, nmbEvents);
181  progress_display* progressIndicator = (not debug) ? new progress_display(nmbEvents, cout, "") : 0;
182  for (long int eventIndex = 0; eventIndex < nmbEvents; ++eventIndex) {
183  if (progressIndicator)
184  ++(*progressIndicator);
185 
186  if ((inTrees[0]->LoadTree(eventIndex) < 0) or (inTrees[1]->LoadTree(eventIndex) < 0))
187  break;
188  inTrees[0]->GetEntry(eventIndex);
189  inTrees[1]->GetEntry(eventIndex);
190 
191  assert( (prodKinMomenta[0]->GetEntriesFast() == nmbProdKinPart[0])
192  and (prodKinMomenta[1]->GetEntriesFast() == nmbProdKinPart[0]));
193  if (debug)
194  printDebug << "event[" << eventIndex << "]: " << nmbProdKinPart[0]
195  << " production kinematics particles:" << endl;
196  prodKinMomenta[2]->Clear();
197  for (int i = 0; i < nmbProdKinPart[0]; ++i)
198  fillDiffLeafs(*(prodKinPartNames[0]), prodKinMomenta, i, absoluteDiff, debug);
199 
200  assert( (decayKinMomenta[0]->GetEntriesFast() == nmbDecayKinPart[0])
201  and (decayKinMomenta[1]->GetEntriesFast() == nmbDecayKinPart[0]));
202  if (debug)
203  printDebug << "event[" << eventIndex << "]: " << nmbDecayKinPart[0]
204  << " decay kinematics particles:" << endl;
205  decayKinMomenta[2]->Clear();
206  for (int i = 0; i < nmbDecayKinPart[0]; ++i)
207  fillDiffLeafs(*(decayKinPartNames[0]), decayKinMomenta, i, absoluteDiff, debug);
208 
209  tree->Fill();
210  if (debug)
211  cout << endl;
212  }
213 
214  tree->OptimizeBaskets(treeCacheSize, 1, "d");
215  tree->Write();
216  outFile->Close();
217  printSucc << "wrote " << nmbEvents << " events to file '" << outFileName << "'" << endl;
218  return true;
219 }
Generated by   doxygen 1.8.1.2