ROOTPWA
plotAmpDiffs.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 // plots differences between amplitude files
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <string>
39 #include <set>
40 
41 #include <boost/progress.hpp>
42 
43 #include "TFile.h"
44 #include "TChain.h"
45 #include "TTreeFormula.h"
46 #include "TCanvas.h"
47 #include "TDirectory.h"
48 #include "TH1.h"
49 #include "TH2.h"
50 
51 #include "reportingUtils.hpp"
52 
53 
54 using namespace std;
55 using namespace boost;
56 using namespace rpwa;
57 
58 
59 bool
60 plotAmpDiffs(const string& inFileNamePattern,
61  const string& outFileName = "ampDiffPlots.root",
62  const long int maxNmbEvents = -1,
63  const string& inTreeName = "ampDiffTree")
64 {
65  // open input file(s)
66  printInfo << "opening input file(s) '" << inFileNamePattern << "'" << endl;
67  TChain inTree(inTreeName.c_str());
68  if (inTree.Add(inFileNamePattern.c_str()) < 1) {
69  printWarn << "no events in input file(s) '" << inFileNamePattern << "'" << endl;
70  return false;
71  }
72  inTree.GetListOfFiles()->ls();
73 
74  // open output file
75  TFile* outFile = 0;
76  if (outFileName != "") {
77  printInfo << "writing difference plots to '" << outFileName << "'" << endl;
78  outFile = TFile::Open(outFileName.c_str(), "RECREATE");
79  }
80 
81  // connect tree leafs
82  TObjString* ampName = 0;
83  UInt_t eventNmb;
84  UInt_t massBinMin, massBinMax; // [MeV/c^2]
85  Double_t valReal[2], valImag[2];
86  Double_t absDiffReal, absDiffImag;
87  Double_t relDiffReal, relDiffImag;
88  inTree.SetBranchAddress("ampName", &ampName);
89  inTree.SetBranchAddress("eventNmb", &eventNmb);
90  inTree.SetBranchAddress("massBinMin", &massBinMin);
91  inTree.SetBranchAddress("massBinMax", &massBinMax);
92  inTree.SetBranchAddress("valReal", valReal);
93  inTree.SetBranchAddress("valImag", valImag);
94  inTree.SetBranchAddress("absDiffReal", &absDiffReal);
95  inTree.SetBranchAddress("absDiffImag", &absDiffImag);
96  inTree.SetBranchAddress("relDiffReal", &relDiffReal);
97  inTree.SetBranchAddress("relDiffImag", &relDiffImag);
98 
99  // book histograms
100  const long int nmbEventsTree = inTree.GetEntries();
101  const long int nmbEvents = ((maxNmbEvents > 0) ? min(maxNmbEvents, nmbEventsTree)
102  : nmbEventsTree);
103  // TH1F* hAmps1Re = new TH1F("amps1Real", "amps1Real;Event Number;#Rgothic[Amp 1]", nmbEvents, -0.5, nmbEvents - 0.5);
104  // TH1F* hAmps1Im = new TH1F("amps1Imag", "amps1Imag;Event Number;#Jgothic[Amp 1]", nmbEvents, -0.5, nmbEvents - 0.5);
105  // TH1F* hAmps2Re = new TH1F("amps2Real", "amps2Real;Event Number;#Rgothic[Amp 2]", nmbEvents, -0.5, nmbEvents - 0.5);
106  // TH1F* hAmps2Im = new TH1F("amps2Imag", "amps2Imag;Event Number;#Jgothic[Amp 2]", nmbEvents, -0.5, nmbEvents - 0.5);
107  TH1F* hAbsDiffRe = new TH1F("absDiffReal", "absDiffReal;#Rgothic[Amp 1 - Amp 2];Count", 1000000, -1e-9, 1e-9);
108  TH1F* hAbsDiffIm = new TH1F("absDiffImag", "absDiffImag;#Jgothic[Amp 1 - Amp 2];Count", 1000000, -1e-8, 1e-8);
109  TH1F* hRelDiffRe = new TH1F("relDiffReal", "relDiffReal;(#Rgothic[Ampl 1] - #Rgothic[Ampl 2]) / #Rgothic[Ampl 1];Count", 1000000, -1e-4, 1e-4);
110  TH1F* hRelDiffIm = new TH1F("relDiffImag", "relDiffImag;(#Jgothic[Ampl 1] - #Jgothic[Ampl 2]) / #Jgothic[Ampl 1];Count", 1000000, -1e-4, 1e-4);
111  // TH2F* hCorrRe = new TH2F("corrReal", "corrReal;#Rgothic[Amp 1];#Rgothic[Amp 2]", 5000, -5, 5, 5000, -5, 5);
112  // TH2F* hCorrIm = new TH2F("corrImag", "corrImag;#Jgothic[Amp 1];#Jgothic[Amp 2]", 5000, -5, 5, 5000, -5, 5);
113 
114  // loop over tree
115  set<string> ampNames;
116  double maxAbsDiff = 0;
117  long int maxAbsDiffIndex = -1;
118  double maxRelDiff = 0;
119  long int maxRelDiffIndex = -1;
120  progress_display progressIndicator(nmbEvents, cout, "");
121  for (long int eventIndex = 0; eventIndex < nmbEvents; ++eventIndex) {
122  ++progressIndicator;
123 
124  if (inTree.LoadTree(eventIndex) < 0)
125  break;
126  inTree.GetEntry(eventIndex);
127 
128  ampNames.insert(ampName->GetString().Data());
129 
130  // hAmps1Re->SetBinContent(eventIndex + 1, valReal[0]);
131  // hAmps1Im->SetBinContent(eventIndex + 1, valImag[0]);
132  // hAmps2Re->SetBinContent(eventIndex + 1, valReal[1]);
133  // hAmps2Im->SetBinContent(eventIndex + 1, valImag[1]);
134  hAbsDiffRe->Fill(absDiffReal);
135  hAbsDiffIm->Fill(absDiffImag);
136  hRelDiffRe->Fill(relDiffReal);
137  hRelDiffIm->Fill(relDiffImag);
138  // hCorrRe->Fill(valReal[0], valReal[1]);
139  // hCorrIm->Fill(valImag[0], valImag[1]);
140  // compute maximum deviations
141  const double absMax = max(fabs(absDiffReal), fabs(absDiffImag));
142  const double relMax = max(fabs(relDiffReal), fabs(relDiffImag));
143  if (absMax > maxAbsDiff) {
144  maxAbsDiff = absMax;
145  maxAbsDiffIndex = eventIndex;
146  }
147  if (relMax > maxRelDiff) {
148  maxRelDiff = relMax;
149  maxRelDiffIndex = eventIndex;
150  }
151  }
152  printInfo << "maximum observed deviations: absolute = " << maxPrecision(maxAbsDiff) << " "
153  << "(event " << maxAbsDiffIndex << "), relative = " << maxPrecision(maxRelDiff) << " "
154  << "(event " << maxRelDiffIndex << ")" << endl;
155 
156  // print events with max deviation
157  inTree.GetEntry(maxAbsDiffIndex);
158  printInfo << "event[" << maxAbsDiffIndex << "]:" << endl
159  << " ampName ....... " << ampName->GetString() << endl
160  << " eventNmb ...... " << eventNmb << endl
161  << " massBinMin .... " << massBinMin << endl
162  << " massBinMax .... " << massBinMax << endl
163  << " valReal[0] .... " << maxPrecision(valReal[0] ) << endl
164  << " valReal[1] .... " << maxPrecision(valReal[1] ) << endl
165  << " valImag[0] .... " << maxPrecision(valImag[0] ) << endl
166  << " valImag[1] .... " << maxPrecision(valImag[1] ) << endl
167  << " absDiffReal ... " << maxPrecision(absDiffReal) << endl
168  << " absDiffImag ... " << maxPrecision(absDiffImag) << endl
169  << " relDiffReal ... " << maxPrecision(relDiffReal) << endl
170  << " relDiffImag ... " << maxPrecision(relDiffImag) << endl;
171  inTree.GetEntry(maxRelDiffIndex);
172  printInfo << "event[" << maxRelDiffIndex << "]:" << endl
173  << " ampName ....... " << ampName->GetString() << endl
174  << " eventNmb ...... " << eventNmb << endl
175  << " massBinMin .... " << massBinMin << endl
176  << " massBinMax .... " << massBinMax << endl
177  << " valReal[0] .... " << maxPrecision(valReal[0] ) << endl
178  << " valReal[1] .... " << maxPrecision(valReal[1] ) << endl
179  << " valImag[0] .... " << maxPrecision(valImag[0] ) << endl
180  << " valImag[1] .... " << maxPrecision(valImag[1] ) << endl
181  << " absDiffReal ... " << maxPrecision(absDiffReal) << endl
182  << " absDiffImag ... " << maxPrecision(absDiffImag) << endl
183  << " relDiffReal ... " << maxPrecision(relDiffReal) << endl
184  << " relDiffImag ... " << maxPrecision(relDiffImag) << endl;
185 
186  const string leafsToDraw[] = {"absDiffReal",
187  "absDiffImag",
188  "relDiffReal",
189  "relDiffImag"};
190  const unsigned int nmbLeafs = sizeof(leafsToDraw) / sizeof(leafsToDraw[0]);
191  TH2F* hLeafs[nmbLeafs];
192  for (unsigned int i = 0; i < nmbLeafs; ++i) {
193  printInfo << "filling histogram for leaf '" << leafsToDraw[i] << "'" << endl;
194  const double min = inTree.GetMinimum(leafsToDraw[i].c_str());
195  const double max = inTree.GetMaximum(leafsToDraw[i].c_str());
196  hLeafs[i] = new TH2F((leafsToDraw[i] + "Wave").c_str(),
197  (leafsToDraw[i] + " vs. Wave;Wave;" + leafsToDraw[i]).c_str(),
198  1, 0, 1, 100000, min - 0.1 * fabs(min), max + 0.1 * fabs(max));
199  hLeafs[i]->SetBit(TH1::kCanRebin);
200  TTreeFormula* leafVal = new TTreeFormula("", leafsToDraw[i].c_str(), &inTree);
201  progress_display progressIndicator(nmbEvents, cout, "");
202  for (long int eventIndex = 0; eventIndex < nmbEvents; ++eventIndex) {
203  ++progressIndicator;
204  if (inTree.LoadTree(eventIndex) < 0)
205  break;
206  inTree.GetEntry(eventIndex);
207  leafVal->UpdateFormulaLeaves();
208  hLeafs[i]->Fill(ampName->GetString().Data(), leafVal->EvalInstance(), 1.);
209  }
210  hLeafs[i]->LabelsDeflate("X");
211  hLeafs[i]->SetLabelSize(0.02, "X");
212  hLeafs[i]->SetDrawOption("COLZ");
213  }
214 
215  if (outFile) {
216  outFile->Write();
217  outFile->Close();
218  delete outFile;
219  printInfo << "wrote difference plots to '" << outFileName << "'" << endl;
220  }
221  return true;
222 }
Generated by   doxygen 1.8.1.2