ROOTPWA
plotSpinTotals.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 // Creates summary plots with intensities of spin totals
27 //
28 //
29 // Author List:
30 // Sebastian Neubert TUM (original author)
31 //
32 //
33 //-----------------------------------------------------------
34 
35 
36 #include <iostream>
37 #include <sstream>
38 #include <vector>
39 
40 #include "TFile.h"
41 #include "TStyle.h"
42 #include "TString.h"
43 #include "TCanvas.h"
44 #include "TROOT.h"
45 #include "TMultiGraph.h"
46 #include "TGraphErrors.h"
47 #include "TAxis.h"
48 #include "TLine.h"
49 #include "TLegend.h"
50 
51 #include "reportingUtils.hpp"
52 #include "fitResult.h"
53 #include "plotSpinTotals.h"
54 
55 
56 using namespace std;
57 using namespace rpwa;
58 
59 
60 vector<pair<string, TVirtualPad*> >
61 plotSpinTotals(const unsigned int nmbTrees, // number of fitResult trees
62  TTree** trees, // array of fitResult trees
63  const int* colors, // array of line and marker colors
64  const double yAxisRangeMax, // if != 0; range of y-axis is not allowed to be larger than this value
65  const bool drawLegend, // if set legend is drawn
66  const string& outFileName,
67  const string& branchName)
68 {
69  vector<pair<string, TVirtualPad*> > wavePads;
70 
71  for (unsigned int i = 0; i < nmbTrees; ++i)
72  if (!trees[i]) {
73  printErr << "null pointer to tree " << i << ". exiting." << endl;
74  return wavePads;
75  }
76 
77  TFile* outFile = NULL;
78  if (outFileName != "")
79  outFile = TFile::Open(outFileName.c_str(), "RECREATE");
80  if (!outFile || outFile->IsZombie()) {
81  printErr << "could not create file '" << outFileName << "'. exiting." << endl;
82  return wavePads;
83  }
84 
85  const unsigned int nmbPadsPerCanvMin = 6; // minimum number of pads each canvas is subdivided into
86  // define set of spin totals
87  const string waves[] = {"logLikelihood",
88  ".*", // total intensity
89  "flat",
90  "0-\\+0\\+",
91  "1\\+\\+0\\+",
92  //"1\\+\\+1\\+",
93  "2-\\+0\\+",
94  "2-\\+1\\+",
95  "2\\+\\+0-",
96  "2\\+\\+1\\+",
97  "1-\\+0-",
98  "1-\\+1\\+",
99  "3\\+\\+0\\+",
100  "4\\+\\+1\\+",
101  "3-\\+1\\+",
102  "3-\\+1-",
103  "3-\\+0-",
104  "0\\+\\+0-"};
105  // const string waves[] = {"logLikelihood",
106  // ".*", // total intensity
107  // "flat",
108  // "0-\\+",
109  // "1\\+\\+",
110  // "1-\\+",
111  // "2\\+\\+",
112  // "2-\\+",
113  // "3\\+\\+",
114  // "4\\+\\+",
115  // "4-\\+",
116  // "0-\\+0\\+",
117  // "1\\+\\+0\\+",
118  // "1\\+\\+1\\+",
119  // "1\\+\\+1-",
120  // "1-\\+0-",
121  // "1-\\+1\\+",
122  // "1-\\+1-",
123  // "2\\+\\+0-",
124  // "2\\+\\+1\\+",
125  // "2\\+\\+1-",
126  // "2-\\+0\\+",
127  // "2-\\+1\\+",
128  // "2-\\+1-",
129  // "3\\+\\+0\\+",
130  // "3\\+\\+1\\+",
131  // "4\\+\\+1\\+",
132  // "4-\\+0\\+",
133  // "4-\\+1\\+"};
134  // const string waves[] = {"logLikelihood",
135  // ".*", // total intensity
136  // "flat",
137  // "1--",
138  // "1\\+-",
139  // "2--",
140  // "3--",
141  // "3\\+-",
142  // "1--0\\+",
143  // "1--1\\+",
144  // "1\\+-1\\+",
145  // "2--1\\+",
146  // "2--2\\+",
147  // "3--0\\+",
148  // "3--1\\+",
149  // "3--2\\+",
150  // "3--3\\+",
151  // "3\\+-1\\+",
152  // "3\\+-2\\+",
153  // "3\\+-3\\+"};
154  const unsigned int nmbWaves = sizeof(waves) / sizeof(string);
155  printInfo << "plotting spin totals for:" << endl;
156  for (unsigned int i = 0; i < nmbWaves; ++i)
157  cout << " " << ((waves[i] != "") ? waves[i] : "total") << endl;
158 
159  // set plot style
160  gStyle->SetOptStat(0);
161  gStyle->SetFillColor(0);
162  gStyle->SetPadColor(0);
163 
164  // calculate optimum canvas subdivision
165  const unsigned int nmbPadsVert = (int)ceil(sqrt(nmbPadsPerCanvMin));
166  const unsigned int nmbPadsHor = (int)ceil((double)nmbPadsPerCanvMin / (double)nmbPadsVert);
167  const unsigned int nmbPadsPerCanv = nmbPadsHor * nmbPadsVert;
168 
169  // plot spin totals
170  wavePads.resize(nmbWaves, pair<string, TVirtualPad*>("", NULL));
171  unsigned int countPad = 0;
172  unsigned int countCanv = 0;
173  TCanvas* canv = 0;
174  string drawOpt;
175  for (unsigned int i = 0; i < nmbWaves; ++i) {
176  const TString waveName = unescapeRegExpSpecialChar((waves[i] != ".*") ? waves[i] : "total");
177  // create new pad, if necessary
178  if (countPad == 0) {
179  stringstream canvName;
180  canvName << "spinTotals" << countCanv;
181  canv = static_cast<TCanvas*>(gROOT->FindObject(canvName.str().c_str()));
182  if (!canv) {
183  canv = new TCanvas(canvName.str().c_str(), "Spin Totals", 10, 10, 1000, 900);
184  canv->Divide(nmbPadsHor, nmbPadsVert);
185  drawOpt = "AP";
186  } else
187  drawOpt = "P SAME";
188  }
189 
190  // build multiGraph
191  TMultiGraph* graph = new TMultiGraph();
192  double maxY = 0;
193  for (unsigned int j = 0; j < nmbTrees; ++j) {
194  // build and run TTree::Draw() expression
195  string drawExpr;
196  if (waves[i] == "logLikelihood") // draw likelihood
197  drawExpr = "-" + branchName + ".logLikelihood() / " + branchName + ".nmbEvents():0:"
198  + branchName + ".massBinCenter()";
199  else
200  drawExpr = branchName + ".intensity(\"" + waves[i] + "\"):"
201  + branchName + ".intensityErr(\"" + waves[i] + "\"):"
202  + branchName + ".massBinCenter()";
203  cout << " running TTree::Draw() expression '" << drawExpr << "' "
204  << "on tree '" << trees[j]->GetName() << "', '" << trees[j]->GetTitle() << "'" << endl;
205  trees[j]->Draw(drawExpr.c_str(), "", "goff");
206 
207  // extract data from TTree::Draw() result and build graph
208  const int nmbBins = trees[j]->GetSelectedRows();
209  vector<double> x(nmbBins), xErr(nmbBins);
210  vector<double> y(nmbBins), yErr(nmbBins);
211  for (int k = 0; k < nmbBins; ++k) {
212  x [k] = trees[j]->GetV3()[k] * 0.001; // convert mass to GeV
213  xErr[k] = 0;
214  y [k] = trees[j]->GetV1()[k];
215  yErr[k] = (waves[i] == "logLikelihood") ? 0 : trees[j]->GetV2()[k];
216  }
217  TGraphErrors* g = new TGraphErrors(nmbBins,
218  &(*(x.begin())), // mass
219  &(*(y.begin())), // intensity
220  &(*(xErr.begin())), // mass error
221  &(*(yErr.begin()))); // intensity error
222 
223  {
224  stringstream graphName;
225  graphName << waveName << "_" << j;
226  g->SetName (graphName.str().c_str());
227  g->SetTitle(graphName.str().c_str());
228  }
229  g->SetMarkerStyle(21);
230  g->SetMarkerSize(0.5);
231  if (colors) {
232  g->SetMarkerColor(colors[j]);
233  g->SetLineColor (colors[j]);
234  }
235  graph->Add(g, "P");
236 
237  // compute maximum for y-axis
238  for (int k = 0; k < nmbBins; ++k) {
239  const double val = y[k] + yErr[k];
240  if (maxY < val)
241  maxY = val;
242  }
243  }
244  cout << " maximum intensity for graph " << graph->GetName() << " is " << maxY << endl;
245 
246  canv->cd(++countPad);
247  graph->SetTitle(waveName);
248  graph->SetName(waveName);
249  if ((yAxisRangeMax > 0) && (maxY > yAxisRangeMax))
250  maxY = yAxisRangeMax;
251  graph->SetMinimum(-maxY * 0.1);
252  graph->SetMaximum( maxY * 1.1);
253  // draw graph
254  graph->Draw(drawOpt.c_str());
255  graph->GetXaxis()->SetTitle("Mass [GeV]");
256  if (waves[i] == "logLikelihood")
257  graph->GetYaxis()->SetTitle("ln(Likelihood)");
258  else
259  graph->GetYaxis()->SetTitle("Intensity");
260  TLine* line = new TLine(graph->GetXaxis()->GetXmin(), 0, graph->GetXaxis()->GetXmax(), 0);
261  line->SetLineStyle(3);
262  graph->GetListOfFunctions()->Add(line);
263  // add legend
264  if (drawLegend && (nmbTrees > 1)) {
265  TLegend* legend = new TLegend(0.65,0.80,0.99,0.99);
266  legend->SetFillColor(10);
267  legend->SetBorderSize(1);
268  legend->SetMargin(0.2);
269  for (unsigned int j = 0; j < nmbTrees; ++j) {
270  TGraph* g = static_cast<TGraph*>(graph->GetListOfGraphs()->At(j));
271  legend->AddEntry(g, trees[j]->GetTitle(), "LPE");
272  }
273  graph->GetListOfFunctions()->Add(legend);
274  }
275 
276  // memorize pad
277  wavePads[i].first = waves[i];
278  wavePads[i].second = gPad;
279 
280  if (outFile)
281  graph->Write();
282 
283  if (countPad == nmbPadsPerCanv) {
284  canv->Update();
285  countPad = 0;
286  ++countCanv;
287  }
288  }
289 
290  if (outFile)
291  outFile->Close();
292 
293  return wavePads;
294 }
Generated by   doxygen 1.8.1.2