ROOTPWA
plotIntensity.C
Go to the documentation of this file.
1 
2 //
3 // Copyright 2009 Sebastian Neubert
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 // $Id$
24 //
25 // Description:
26 // Draws intensity graph for single wave from tree.
27 //
28 //
29 // Author List:
30 // Sebastian Neubert TUM (original author)
31 //
32 //
33 //-----------------------------------------------------------
34 
35 
36 #include <iostream>
37 #include <sstream>
38 
39 #include "TGraphErrors.h"
40 #include "TAxis.h"
41 #include "TLine.h"
42 #include "TPad.h"
43 #include "TLegend.h"
44 #include "TList.h"
45 
46 #include "reportingUtils.hpp"
47 #include "plotIntensity.h"
48 
49 
50 using namespace std;
51 using namespace rpwa;
52 
53 // signature with wave name
54 TMultiGraph*
55 plotIntensity(const unsigned int nmbTrees, // number of fitResult trees
56  TTree** trees, // array of fitResult trees
57  const std::string& waveName, // wave index
58  const bool saveEps, // if set, EPS file with name wave ID
59  const int* graphColors, // array of colors for graph line and marker
60  const double* graphScales, // array of scales for graphgroups
61  const bool drawLegend, // if set legend is drawn
62  const string& graphTitle, // name and title of graph (default is wave IDs)
63  const char* drawOption, // draw option for graph
64  const double normalization, // scale factor for intensities
65  const double yAxisRangeMax, // if != 0; range of y-axis is limited to this value
66  const string& selectExpr, // TTree::Draw() selection expression
67  const string& branchName) // fitResult branch name
68 {
69 
70  for (unsigned int i = 0; i < nmbTrees; ++i)
71  if (!trees[i]) {
72  printErr << "null pointer to tree[" << i << "]. aborting." << endl;
73  return 0;
74  }
75  printInfo << "plotting wave intensity for wave '" << waveName;
76  if (selectExpr != "")
77  cout << " using selection criterion '" << selectExpr << "'";
78  cout << endl;
79 
80  // create multiGraph
81  TMultiGraph* graph = new TMultiGraph();
82  {
83  if (graphTitle != "") {
84  graph->SetName (graphTitle.c_str());
85  graph->SetTitle(graphTitle.c_str());
86  } else {
87  stringstream suffix;
88  //suffix << " [" << "no waveIndex" << "]";
89  graph->SetName (waveName.c_str());
90  graph->SetTitle((waveName + suffix.str()).c_str());
91  }
92  }
93 
94  // fill multiGraph
95  double maxY = 0;
96  for (unsigned int i = 0; i < nmbTrees; ++i) {
97 
98  // get wave index for this tree
99  fitResult* massBin = new fitResult();
100  trees[i]->SetBranchAddress(branchName.c_str(), &massBin);
101  vector<int> waveIndexThisTree; // vector of wave indices per bin
102  bool found(false);
103  for (int imassBin = 0; imassBin < trees[i]->GetEntries(); imassBin++){
104  trees[i]->GetEntry(imassBin);
105  waveIndexThisTree.push_back(massBin->waveIndex(waveName));
106  if (waveIndexThisTree[imassBin] >= 0){
107  found = true;
108  }
109  }
110  if (!found) {
111  printInfo << "cannot find wave '" << waveName << "' in tree '" << trees[i]->GetTitle() << "'. "
112  << "skipping." << endl;
113  continue;
114  }
115 
116  // running bin by bin since wave name does necessarily correspond to a fixed wave number
117 
118  // build and run TTree::Draw() expression
119  stringstream drawExpr;
120  drawExpr << branchName << ".intensity(\"" << waveName << "\"):"
121  << branchName << ".intensityErr(\"" << waveName << "\"):"
122  << branchName << ".massBinCenter() >> h" << waveName << "_" << i;
123  cout << " running TTree::Draw() expression '" << drawExpr.str() << "' "
124  << "on tree '" << trees[i]->GetName() << "', '" << trees[i]->GetTitle() << "'" << endl;
125  trees[i]->Draw(drawExpr.str().c_str(), selectExpr.c_str(), "goff");
126 
127  // extract data from TTree::Draw() result and build graph
128  const int nmbBins = trees[i]->GetSelectedRows();
129  vector<double> x, xErr;
130  vector<double> y, yErr;
131 
132  double scale = normalization;
133  if(graphScales) {
134  scale = scale * graphScales[i];
135  }
136 
137  for (int j = 0; j < nmbBins; ++j) {
138  y .push_back(trees[i]->GetV1()[j] * scale); // scale intensities
139  // if (y[y.size()-1] == 0){ // remove 0 entries
140  // y.pop_back();
141  // continue;
142  // }
143  x .push_back(trees[i]->GetV3()[j] * 0.001); // convert mass to GeV
144  xErr.push_back(0);
145  yErr.push_back(trees[i]->GetV2()[j] * scale); // scale intensity errors
146  }
147  TGraphErrors* g = new TGraphErrors(x.size(),
148  &(*(x.begin())), // mass
149  &(*(y.begin())), // intensity
150  &(*(xErr.begin())), // mass error
151  &(*(yErr.begin()))); // intensity error
152 
153  // beautify graph
154  stringstream graphName;
155  graphName << graph->GetName() << "_" << i;
156  g->SetName (graphName.str().c_str());
157  g->SetTitle(graphName.str().c_str());
158  g->SetMarkerStyle(21);
159  g->SetMarkerSize(0.5);
160  if (graphColors) {
161  g->SetMarkerColor(graphColors[i]);
162  g->SetLineColor (graphColors[i]);
163  }
164  graph->Add(g);
165 
166  // compute maximum for y-axis
167  for (unsigned int j = 0; j < x.size(); ++j) {
168  const double val = y[j] + yErr[j];
169  if (maxY < val)
170  maxY = val;
171  }
172  }
173  cout << " maximum intensity for graph " << graph->GetName() << " is " << maxY << endl;
174 
175  if ((yAxisRangeMax > 0) && (maxY > yAxisRangeMax))
176  maxY = yAxisRangeMax;
177  graph->SetMinimum(-maxY * 0.1);
178  graph->SetMaximum( maxY * 1.1);
179  // draw graph
180  graph->Draw(drawOption);
181  graph->GetXaxis()->SetTitle("Mass [GeV]");
182  graph->GetYaxis()->SetTitle("Intensity");
183  TLine line;
184  line.SetLineStyle(3);
185  line.DrawLine(graph->GetXaxis()->GetXmin(), 0, graph->GetXaxis()->GetXmax(), 0);
186  gPad->Update();
187 
188  // add legend
189  if (drawLegend && (nmbTrees > 1)) {
190  TLegend* legend = new TLegend(0.65,0.80,0.99,0.99);
191  legend->SetFillColor(10);
192  legend->SetBorderSize(1);
193  legend->SetMargin(0.2);
194  for (unsigned int i = 0; i < nmbTrees; ++i) {
195  TGraph* g = static_cast<TGraph*>(graph->GetListOfGraphs()->At(i));
196  legend->AddEntry(g, trees[i]->GetTitle(), "LPE");
197  }
198  legend->Draw();
199  }
200 
201  // create EPS file
202  if (saveEps)
203  gPad->SaveAs(((string)graph->GetName() + ".eps").c_str());
204 
205  return graph;
206 }
207 
208 // signature with wave index
209 TMultiGraph*
210 plotIntensity(const unsigned int nmbTrees, // number of fitResult trees
211  TTree** trees, // array of fitResult trees
212  const int waveIndex, // wave index
213  const bool saveEps, // if set, EPS file with name wave ID
214  const int* graphColors, // array of colors for graph line and marker
215  const double* graphScales, // array of scales for graphgroups
216  const bool drawLegend, // if set legend is drawn
217  const string& graphTitle, // name and title of graph (default is wave IDs)
218  const char* drawOption, // draw option for graph
219  const double normalization, // scale factor for intensities
220  const double yAxisRangeMax, // if != 0; range of y-axis is limited to this value
221  const string& selectExpr, // TTree::Draw() selection expression
222  const string& branchName) // fitResult branch name
223 {
224  for (unsigned int i = 0; i < nmbTrees; ++i)
225  if (!trees[i]) {
226  printErr << "null pointer to tree[" << i << "]. aborting." << endl;
227  return 0;
228  }
229  string waveName;
230  {
231  // get wave name (assume same wave set in all bins)
232  fitResult* massBin = new fitResult();
233  trees[0]->SetBranchAddress(branchName.c_str(), &massBin);
234  trees[0]->GetEntry(0);
235  waveName = massBin->waveName(waveIndex).Data();
236  }
237  printInfo << "plotting wave intensity for wave '" << waveName << "' [" << waveIndex << "]";
238  if (selectExpr != "")
239  cout << " using selection criterion '" << selectExpr << "'";
240  cout << endl;
241 
242  // create multiGraph
243  TMultiGraph* graph = new TMultiGraph();
244  {
245  if (graphTitle != "") {
246  graph->SetName (graphTitle.c_str());
247  graph->SetTitle(graphTitle.c_str());
248  } else {
249  stringstream suffix;
250  suffix << " [" << waveIndex << "]";
251  graph->SetName (waveName.c_str());
252  graph->SetTitle((waveName + suffix.str()).c_str());
253  }
254  }
255 
256  // fill multiGraph
257  double maxY = 0;
258  for (unsigned int i = 0; i < nmbTrees; ++i) {
259 
260  // get wave index for this tree (assume same wave set in all bins)
261  fitResult* massBin = new fitResult();
262  trees[i]->SetBranchAddress(branchName.c_str(), &massBin);
263  trees[i]->GetEntry(0);
264  const int waveIndexThisTree = massBin->waveIndex(waveName);
265  if (waveIndexThisTree < 0) {
266  printInfo << "cannot find wave '" << waveName << "' in tree '" << trees[i]->GetTitle() << "'. "
267  << "skipping." << endl;
268  continue;
269  }
270 
271  // build and run TTree::Draw() expression
272  stringstream drawExpr;
273  drawExpr << branchName << ".intensity(" << waveIndexThisTree << "):"
274  << branchName << ".intensityErr(" << waveIndexThisTree << "):"
275  << branchName << ".massBinCenter() >> h" << waveName << "_" << i;
276  cout << " running TTree::Draw() expression '" << drawExpr.str() << "' "
277  << "on tree '" << trees[i]->GetName() << "', '" << trees[i]->GetTitle() << "'" << endl;
278  trees[i]->Draw(drawExpr.str().c_str(), selectExpr.c_str(), "goff");
279 
280  // extract data from TTree::Draw() result and build graph
281  const int nmbBins = trees[i]->GetSelectedRows();
282  vector<double> x(nmbBins), xErr(nmbBins);
283  vector<double> y(nmbBins), yErr(nmbBins);
284 
285  double scale = normalization;
286  if (graphScales) {
287  scale = scale * graphScales[i];
288  }
289 
290  for (int j = 0; j < nmbBins; ++j) {
291  x [j] = trees[i]->GetV3()[j] * 0.001; // convert mass to GeV
292  xErr[j] = 0;
293  y [j] = trees[i]->GetV1()[j] * scale; // scale intensities
294  yErr[j] = trees[i]->GetV2()[j] * scale; // scale intensity errors
295  }
296  TGraphErrors* g = new TGraphErrors(nmbBins,
297  &(*(x.begin())), // mass
298  &(*(y.begin())), // intensity
299  &(*(xErr.begin())), // mass error
300  &(*(yErr.begin()))); // intensity error
301 
302  // beautify graph
303  stringstream graphName;
304  graphName << graph->GetName() << "_" << i;
305  g->SetName (graphName.str().c_str());
306  g->SetTitle(graphName.str().c_str());
307  g->SetMarkerStyle(21);
308  g->SetMarkerSize(0.5);
309  if (graphColors) {
310  g->SetMarkerColor(graphColors[i]);
311  g->SetLineColor (graphColors[i]);
312  }
313  graph->Add(g);
314 
315  // compute maximum for y-axis
316  for (int j = 0; j < nmbBins; ++j) {
317  const double val = y[j] + yErr[j];
318  if (maxY < val)
319  maxY = val;
320  }
321  }
322  cout << " maximum intensity for graph " << graph->GetName() << " is " << maxY << endl;
323 
324  if ((yAxisRangeMax > 0) && (maxY > yAxisRangeMax))
325  maxY = yAxisRangeMax;
326  graph->SetMinimum(-maxY * 0.1);
327  graph->SetMaximum( maxY * 1.1);
328  // draw graph
329  graph->Draw(drawOption);
330  graph->GetXaxis()->SetTitle("Mass [GeV]");
331  graph->GetYaxis()->SetTitle("Intensity");
332  TLine line;
333  line.SetLineStyle(3);
334  line.DrawLine(graph->GetXaxis()->GetXmin(), 0, graph->GetXaxis()->GetXmax(), 0);
335  gPad->Update();
336 
337  // add legend
338  if (drawLegend && (nmbTrees > 1)) {
339  TLegend* legend = new TLegend(0.65,0.80,0.99,0.99);
340  legend->SetFillColor(10);
341  legend->SetBorderSize(1);
342  legend->SetMargin(0.2);
343  for (unsigned int i = 0; i < nmbTrees; ++i) {
344  TGraph* g = static_cast<TGraph*>(graph->GetListOfGraphs()->At(i));
345  legend->AddEntry(g, trees[i]->GetTitle(), "LPE");
346  }
347  legend->Draw();
348  }
349 
350  // create EPS file
351  if (saveEps)
352  gPad->SaveAs(((string)graph->GetName() + ".eps").c_str());
353 
354  return graph;
355 }
Generated by   doxygen 1.8.1.2