ROOTPWA
plotAllIntensities.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 graphs for all waves.
27 //
28 //
29 // Author List:
30 // Sebastian Neubert TUM (original author)
31 //
32 //
33 //-----------------------------------------------------------
34 
35 
36 #include <iostream>
37 #include <sstream>
38 #include <map>
39 #include <algorithm>
40 #include <cmath>
41 
42 #include "TCanvas.h"
43 #include "TROOT.h"
44 #include "TFile.h"
45 #include "TGraph.h"
46 #include "TMultiGraph.h"
47 #include "TLatex.h"
48 #include "TList.h"
49 #include "TPostScript.h"
50 #include "TSystem.h"
51 
52 #include "fitResult.h"
53 #include "plotIntensity.h"
54 #include "plotAllIntensities.h"
55 
56 
57 using namespace std;
58 using namespace rpwa;
59 
60 
61 string
62 buildCanvName(const string& jpc,
63  const int index)
64 {
65  string name;
66  {
67  stringstream n;
68  n << jpc << index;
69  name = n.str();
70  }
71  // remove all spaces
72  string::size_type pos = name.find(" ");
73  while (pos != string::npos) {
74  name.replace(pos, 1, "");
75  pos = name.find(" ", pos + 1);
76  }
77  return name;
78 }
79 
80 
81 // predicate for sort
82 bool
83 compareIntensities(const pair<string, double>& a,
84  const pair<string, double>& b)
85 {
86  return a.second > b.second;
87 }
88 
89 vector<pair<string, TVirtualPad*> >
90 plotAllIntensities(const unsigned int nmbTrees, // number of fitResult trees
91  TTree** trees, // array of fitResult trees
92  const bool createPsFile, // if true, plots are written to waves.ps
93  const string& outPath, // path for output files
94  const int* graphColors, // array of colors for graph line and marker
95  const double* graphScales, // array of colors for graph line and marker
96  const bool drawLegend, // if set legend is drawn
97  const double yAxisRangeMax, // if != 0; range of y-axis is limited to this value
98  const string& branchName)
99 {
100  const double intensityThr = 0; // threshold for total intensity in mass bin
101  const int nmbPadsPerCanvMin = 4; // minimum number of pads each canvas is subdivided into
102  vector<pair<string, TVirtualPad*> > wavePads; // return value
103 
104  fitResult* massBin = new fitResult();
105  map<string, int> wavelist; // map with counts of all available waves
106  double totIntensity = 0; // total intensity of all waves
107  int nmbBinsAboveThr = 0;
108  map<string, double> totIntensities; // total intensities per wave name
109  int nbinelements(0);
110 
111  for (unsigned int i = 0; i < nmbTrees; ++i)
112  if (!trees[i]) {
113  printErr << "null pointer to tree " << i << ". exiting." << endl;
114  return wavePads;
115  } else {
116  // check whether all mass bins have the same wave set and store the maximum count
117  trees[i]->SetBranchAddress(branchName.c_str(), &massBin);
118  for (int imassbin = 0; imassbin < trees[i]->GetEntries(); imassbin++){
119  trees[i]->GetEntry(imassbin);
120  nbinelements++;
121  const double binIntensity = massBin->intensity();
122  totIntensity += binIntensity;
123  for (unsigned iwave = 0; iwave < massBin->waveNames().size(); iwave++){
124  wavelist[massBin->waveNames()[iwave]]++;
125  // found a new wave when the specific counter is 1
126  if (wavelist[massBin->waveNames()[iwave]] == 1){
127 
128  }
129  // calculate the total wave intensities
130  // warning, several fit results per bin and wave name lead to a wrong calculation
131  if (binIntensity > intensityThr) {
132  totIntensities[massBin->waveNames()[iwave]]+=massBin->intensity(iwave);
133  if (iwave == 0)
134  ++nmbBinsAboveThr;
135  }
136  }
137  }
138  }
139  //const int nmbWaves = wavelist.size();
140  printInfo << "drawing wave intensities from tree '" << trees[0]->GetName() << "' for "
141  << wavelist.size() << " waves in " << nbinelements << " mass bins." << endl;
142 
143  printInfo << "sorting waves according to their intensity" << endl;
144  vector<pair<string, double> > waveIntensities(totIntensities.size(), pair<string, double>("", 0));
145  {
146  int i(0);
147  for (map<string,double>::iterator it_intensity = totIntensities.begin(); it_intensity != totIntensities.end(); it_intensity++){
148  const double relIntensity = (*it_intensity).second / totIntensity;
149  const string waveName = (*it_intensity).first;
150  waveIntensities[i] = make_pair(waveName, relIntensity);
151  i++;
152  }
153  }
154  sort(waveIntensities.begin(), waveIntensities.end(), compareIntensities); // map is not suitable for sorting according to second term
155 
156  printInfo << "grouping waves w.r.t. to their JPC ..." << endl;
157  map<string, int> nmbWavesPerJpc; // number of waves for each JPC
158  {
159  int i(0);
160  for (map<string,double>::iterator it_intensity = totIntensities.begin(); it_intensity != totIntensities.end(); it_intensity++){
161  const string waveName = (*it_intensity).first;
162  const string jpc = (waveName != "flat") ? string(waveName, 2, 3) : waveName;
163  cout << " wave [" << i << "] of " << totIntensities.size() << ": " << waveName << ", JPC = " << jpc << endl;
164  ++nmbWavesPerJpc[jpc];
165  i++;
166  }
167  }
168  cout << " ... finished" << endl;
169 
170  printInfo << "creating canvases for all JPCs ..." << endl;
171  const int nmbPadsHor = (int)ceil(sqrt(nmbPadsPerCanvMin));
172  const int nmbPadsVert = (int)ceil((double)nmbPadsPerCanvMin / (double)nmbPadsHor);
173  const int nmbPadsPerCanv = nmbPadsHor * nmbPadsVert;
174  map<string, TCanvas*> canvases;
175  for (map<string, int>::const_iterator i = nmbWavesPerJpc.begin(); i != nmbWavesPerJpc.end(); ++i) {
176  const string jpc = i->first;
177  const int nmbCanv = (int)ceil((double)i->second / (double)nmbPadsPerCanv);
178  for (int j = 0; j < nmbCanv; ++j) {
179  const string name = buildCanvName(jpc, j);
180  if (canvases[name] == NULL) {
181  canvases[name] = new TCanvas(name.c_str(), name.c_str(), 10, 10, 1000, 800);
182  canvases[name]->Divide(nmbPadsHor, nmbPadsVert);
183  }
184  }
185  }
186 
187  printInfo << "creating output ROOT file ..." << endl;
188  const string outFileName = outPath + "waveIntensities.root";
189  TFile* outFile = TFile::Open(outFileName.c_str(), "RECREATE");
190  if (!outFile) {
191  printErr << "could not create output file '" << outFileName << "'. exiting." << endl;
192  return wavePads;
193  }
194  gROOT->cd();
195 
196  printInfo << "plotting waves ordered by decreasing intensity ..." << endl;
197  wavePads.resize(waveIntensities.size(), pair<string, TVirtualPad*>("", NULL));
198  map<string, int> canvJpcCounter;
199  for (unsigned int i = 0; i < waveIntensities.size(); ++i) {
200  // get canvas
201  const string waveName = waveIntensities[i].first;
202  const string jpc = (waveName != "flat") ? string(waveName, 2, 3) : waveName;
203  const int canvCount = canvJpcCounter[jpc];
204  const string canvName = buildCanvName(jpc, (int)floor((double)canvCount / (double)nmbPadsPerCanv));
205  const int padIndex = (canvCount % nmbPadsPerCanv) + 1;
206  TCanvas* canv = canvases[canvName];
207  canv->cd(padIndex);
208  cout << endl
209  << " Selecting canvas '" << canvName << "'." << endl
210  << " wave " << canvCount + 1 << " of " << nmbWavesPerJpc[jpc]
211  << " for JPC = " << jpc << ", intensity rank = " << i + 1 << endl;
212  ++canvJpcCounter[jpc];
213 
214  // draw intensity graph
215  TMultiGraph* graph = plotIntensity(nmbTrees, trees, waveName, false, graphColors, graphScales, drawLegend,
216  "", "AP", 1, yAxisRangeMax, "", branchName);
217  if (!graph)
218  continue;
219 
220  // write graph to file
221  outFile->cd();
222  graph->Write();
223  gROOT->cd();
224 
225  // draw additional info
226  const double intensity = floor(waveIntensities[i].second * 1000 + 0.5) / 10;
227  stringstream label;
228  label << "I = " << intensity << "% (" << i + 1 << ")";
229  TLatex* text = new TLatex(0.15, 0.85, label.str().c_str());
230  text->SetNDC(true);
231  text->Draw();
232 
233  // memorize pad
234  wavePads[i].first = waveName;
235  wavePads[i].second = gPad;
236  }
237 
238  // also write TH3s created by plotIntensity() to file
239  TList* histList = gDirectory->GetList();
240  for (unsigned int i = 0; i < nmbTrees; ++i)
241  histList->Remove(trees[i]);
242  printInfo << "writing the following " << histList->GetEntries() <<" objects to '" << outFile->GetName() << "'" << endl;
243  histList->Print();
244  outFile->cd();
245  histList->Write();
246 
247  // cleanup
248  outFile->Close();
249  if (outFile)
250  delete outFile;
251  outFile = 0;
252  gROOT->cd();
253 
254  if (createPsFile) {
255  const string psFileName = outPath + "waveIntensities.ps";
256  TCanvas dummyCanv("dummy", "dummy");
257  dummyCanv.Print((psFileName + "[").c_str());
258  for (map<string, TCanvas*>::iterator i = canvases.begin(); i != canvases.end(); ++i) {
259  i->second->Print(psFileName.c_str());
260  delete i->second;
261  i->second = 0;
262  }
263  dummyCanv.Print((psFileName + "]").c_str());
264  gSystem->Exec(("gv " + psFileName).c_str());
265  }
266 
267  return wavePads;
268 }
269 
270 
271 /*
272 vector<pair<string, TVirtualPad*> >
273 plotAllIntensities(const unsigned int nmbTrees, // number of fitResult trees
274  TTree** trees, // array of fitResult trees
275  const bool createPsFile, // if true, plots are written to waves.ps
276  const string& outPath, // path for output files
277  const int* graphColors, // array of colors for graph line and marker
278  const bool drawLegend, // if set legend is drawn
279  const double yAxisRangeMax, // if != 0; range of y-axis is limited to this value
280  const string& branchName)
281 {
282  const double intensityThr = 0; // threshold for total intensity in mass bin
283  const int nmbPadsPerCanvMin = 4; // minimum number of pads each canvas is subdivided into
284  vector<pair<string, TVirtualPad*> > wavePads; // return value
285 
286  for (unsigned int i = 0; i < nmbTrees; ++i)
287  if (!trees[i]) {
288  printErr << "null pointer to tree " << i << ". exiting." << endl;
289  return wavePads;
290  }
291  // assume that all mass trees have same wave set
292  fitResult* massBin = new fitResult();
293  trees[0]->SetBranchAddress(branchName.c_str(), &massBin);
294  const int nmbMassBins = trees[0]->GetEntries();
295  trees[0]->GetEntry(0);
296  const int nmbWaves = massBin->nmbWaves(); // assumes that number of waves is the same for all bins
297  printInfo << "drawing wave intensities from tree '" << trees[0]->GetName() << "' for "
298  << nmbWaves << " waves in " << nmbMassBins << " mass bins." << endl;
299 
300  printInfo << "calculating total wave intensities" << endl;
301  int nmbBinsAboveThr = 0;
302  double totIntensity = 0; // total intensity of all waves
303  vector<double> totIntensities(nmbWaves, 0); // total intensity of individual waves
304  for (int i = 0; i < nmbMassBins; ++i) {
305  trees[0]->GetEntry(i);
306  const double binIntensity = massBin->intensity();
307  totIntensity += binIntensity;
308  if (binIntensity > intensityThr) {
309  for (int j = 0; j < nmbWaves; ++j)
310  totIntensities[j] += massBin->intensity(j);
311  ++nmbBinsAboveThr;
312  }
313  }
314 
315  printInfo << "sorting waves according to their intensity" << endl;
316  vector<pair<string, double> > waveIntensities(nmbWaves, pair<string, double>("", 0));
317  for (int i = 0; i < nmbWaves; ++i) {
318  const double relIntensity = totIntensities[i] / totIntensity;
319  const string waveName = massBin->waveName(i).Data();
320  waveIntensities[i] = make_pair(waveName, relIntensity);
321  }
322  sort(waveIntensities.begin(), waveIntensities.end(), compareIntensities);
323 
324  printInfo << "grouping waves w.r.t. to their JPC ..." << endl;
325  map<string, int> nmbWavesPerJpc; // number of waves for each JPC
326  for (int i = 0; i < nmbWaves; ++i) {
327  const string waveName = massBin->waveName(i).Data();
328  const string jpc = (waveName != "flat") ? string(waveName, 2, 3) : waveName;
329  cout << " wave [" << i << "] of " << nmbWaves << ": " << waveName << ", JPC = " << jpc << endl;
330  ++nmbWavesPerJpc[jpc];
331  }
332  cout << " ... finished" << endl;
333 
334  printInfo << "creating canvases for all JPCs ..." << endl;
335  const int nmbPadsHor = (int)ceil(sqrt(nmbPadsPerCanvMin));
336  const int nmbPadsVert = (int)ceil((double)nmbPadsPerCanvMin / (double)nmbPadsHor);
337  const int nmbPadsPerCanv = nmbPadsHor * nmbPadsVert;
338  map<string, TCanvas*> canvases;
339  for (map<string, int>::const_iterator i = nmbWavesPerJpc.begin(); i != nmbWavesPerJpc.end(); ++i) {
340  const string jpc = i->first;
341  const int nmbCanv = (int)ceil((double)i->second / (double)nmbPadsPerCanv);
342  for (int j = 0; j < nmbCanv; ++j) {
343  const string name = buildCanvName(jpc, j);
344  if (canvases[name] == NULL) {
345  canvases[name] = new TCanvas(name.c_str(), name.c_str(), 10, 10, 1000, 800);
346  canvases[name]->Divide(nmbPadsHor, nmbPadsVert);
347  }
348  }
349  }
350 
351  printInfo << "creating output ROOT file ..." << endl;
352  const string outFileName = outPath + "waveIntensities.root";
353  TFile* outFile = TFile::Open(outFileName.c_str(), "RECREATE");
354  if (!outFile) {
355  printErr << "could not create output file '" << outFileName << "'. exiting." << endl;
356  return wavePads;
357  }
358  gROOT->cd();
359 
360  printInfo << "plotting waves ordered by decreasing intensity ..." << endl;
361  wavePads.resize(nmbWaves, pair<string, TVirtualPad*>("", NULL));
362  map<string, int> canvJpcCounter;
363  for (int i = 0; i < nmbWaves; ++i) {
364  // get canvas
365  const string waveName = waveIntensities[i].first;
366  const string jpc = (waveName != "flat") ? string(waveName, 2, 3) : waveName;
367  const int canvCount = canvJpcCounter[jpc];
368  const string canvName = buildCanvName(jpc, (int)floor((double)canvCount / (double)nmbPadsPerCanv));
369  const int padIndex = (canvCount % nmbPadsPerCanv) + 1;
370  TCanvas* canv = canvases[canvName];
371  canv->cd(padIndex);
372  cout << endl
373  << " Selecting canvas '" << canvName << "'." << endl
374  << " wave " << canvCount + 1 << " of " << nmbWavesPerJpc[jpc]
375  << " for JPC = " << jpc << ", intensity rank = " << i + 1 << endl;
376  ++canvJpcCounter[jpc];
377 
378  // draw intensity graph
379  TMultiGraph* graph = plotIntensity(nmbTrees, trees, waveName, false, graphColors, drawLegend,
380  "", "AP", 1, yAxisRangeMax, "", branchName);
381  if (!graph)
382  continue;
383 
384  // write graph to file
385  outFile->cd();
386  graph->Write();
387  gROOT->cd();
388 
389  // draw additional info
390  const double intensity = floor(waveIntensities[i].second * 1000 + 0.5) / 10;
391  stringstream label;
392  label << "I = " << intensity << "% (" << i + 1 << ")";
393  TLatex* text = new TLatex(0.15, 0.85, label.str().c_str());
394  text->SetNDC(true);
395  text->Draw();
396 
397  // memorize pad
398  wavePads[i].first = waveName;
399  wavePads[i].second = gPad;
400  }
401 
402  // also write TH3s created by plotIntensity() to file
403  TList* histList = gDirectory->GetList();
404  for (unsigned int i = 0; i < nmbTrees; ++i)
405  histList->Remove(trees[i]);
406  printInfo << "writing the following " << histList->GetEntries() <<" objects to '" << outFile->GetName() << "'" << endl;
407  histList->Print();
408  outFile->cd();
409  histList->Write();
410 
411  // cleanup
412  outFile->Close();
413  if (outFile)
414  delete outFile;
415  outFile = 0;
416  gROOT->cd();
417 
418  if (createPsFile) {
419  const string psFileName = outPath + "waveIntensities.ps";
420  TCanvas dummyCanv("dummy", "dummy");
421  dummyCanv.Print((psFileName + "[").c_str());
422  for (map<string, TCanvas*>::iterator i = canvases.begin(); i != canvases.end(); ++i) {
423  i->second->Print(psFileName.c_str());
424  delete i->second;
425  i->second = 0;
426  }
427  dummyCanv.Print((psFileName + "]").c_str());
428  gSystem->Exec(("gv " + psFileName).c_str());
429  }
430 
431  return wavePads;
432 }
433 */
Generated by   doxygen 1.8.1.2