ROOTPWA
compareWithDima.C
Go to the documentation of this file.
1 //
2 // plots result from Dima's fit on top of ROOTPWA histograms
3 //
4 
5 
6 #include <string>
7 #include <iostream>
8 #include <map>
9 #include <algorithm>
10 #include <vector>
11 #include <utility>
12 
13 #include "TFile.h"
14 #include "TMultiGraph.h"
15 #include "TGraphErrors.h"
16 #include "TRegexp.h"
17 #include "TPRegexp.h"
18 #include "TH1.h"
19 #include "TLegend.h"
20 #include "TKey.h"
21 #include "TObject.h"
22 #include "TString.h"
23 #include "TCanvas.h"
24 #include "TPad.h"
25 #include "TPostScript.h"
26 #include "TSystem.h"
27 
28 #include "reportingUtils.hpp"
29 #include "fitResult.h"
30 #include "plotAllIntensities.h"
31 #include "plotIntensity.h"
32 #include "plotSpinTotals.h"
33 #include "plotPhase.h"
34 
35 
36 using namespace std;
37 using namespace rpwa;
38 
39 
40 // reads Dima's histograms from file into map (histogram title, histogram pointer)
41 map<string, TH1*>
42 readDimaHists(TFile* dimaFile,
43  TPRegexp histNamePattern) // histogram name pattern
44 {
45  map<string, TH1*> dimaHists;
46  TIterator* keys = dimaFile->GetListOfKeys()->MakeIterator();
47  while (TKey* k = static_cast<TKey*>(keys->Next())) {
48  if (!k) {
49  printWarn << "NULL pointer to TKey. skipping." << endl;
50  continue;
51  }
52  TObject* o = k->ReadObj();
53  if (!o) {
54  printWarn << "cannot read object from TKey '" << k->GetName() << "'. skipping." << endl;
55  continue;
56  }
57  const string type = o->IsA()->GetName();
58  if (type.find("TH1F") != string::npos) {
59  TH1* h = static_cast<TH1*>(o);
60  const TString hName = h->GetName();
61  if (hName.Contains(histNamePattern))
62  dimaHists[h->GetTitle()] = h;
63  }
64  }
65  return dimaHists;
66 }
67 
68 
69 // constructs wave name in Dima's convention from wave name in ROOTPWA convention
70 string translateWaveName(const string& waveName)
71 {
72  TString dimaName = waveName;
73 
74  if (dimaName == "flat")
75  return "FLAT";
76 
77  // add braces and space
78  dimaName.Insert(2, "(");
79  dimaName.Insert(6, ")");
80  dimaName.Insert(9, " ");
81 
82  // translate isobars
83  {
84  const unsigned int nmbDict = 5;
85  const string dictionary[nmbDict][2] = {//{"sigma", "eps1"},
86  {"sigma", "f0(1400)"},
87  //{"f0980", "f0"},
88  {"f0980", "f0(980)"},
89  {"rho770", "rho"},
90  {"f21270", "f2"},
91  {"rho31690", "rho3"}};
92  for (unsigned int i = 0; i < nmbDict; ++i)
93  dimaName.ReplaceAll(dictionary[i][0].c_str(), dictionary[i][1].c_str());
94  }
95 
96  // translate L
97  {
98  map<char, char> dictionary;
99  dictionary['0'] = 'S';
100  dictionary['1'] = 'P';
101  dictionary['2'] = 'D';
102  dictionary['3'] = 'F';
103  dictionary['4'] = 'G';
104  dictionary['5'] = 'H';
105  dictionary['6'] = 'I';
106  dictionary['7'] = 'K';
107  dictionary['8'] = 'L';
108 
109  int pos = dimaName.First("_");
110  char L = dimaName(pos + 1);
111 
112  // cut off last part
113  dimaName = dimaName(0, pos);
114  dimaName += " pi ";
115  dimaName += dictionary[L];
116  }
117 
118  return dimaName.Data();
119 }
120 
121 
122 void
123 compareIntensityWithDima(TTree* tree, TFile* dimaFile, std::string waveName) {
124  printInfo << "comparing intensities with Dima's result for wave " << waveName << endl;
125 
126  // get Dima's intensity histograms from file
127  map<string, TH1*> dimaIntensityHists; // intensity histograms w/o acceptance correction (see histgroups.txt)
128  {
129  map<string, TH1*> dimaHists = readDimaHists(dimaFile, TString("^h10\\d{2}$"));
130  // transform map with histogram titles as keys to map with wave names as keys
131  for (map<string, TH1*>::iterator i = dimaHists.begin(); i != dimaHists.end(); ++i) {
132  cout << " found intensity histogram '"<< i->first << "' " << flush;
133  TString waveName = i->first;
134  waveName = waveName.Strip(TString::kTrailing, ' ');
135  dimaIntensityHists[waveName.Data()] = i->second;
136  cout << "for wave '" << waveName << "'" << endl;
137  }
138  }
139 
140  TCanvas *pad = new TCanvas(waveName.c_str());
141  pad->cd();
142  TMultiGraph* rootgraph = plotIntensity(tree, waveName);
143  rootgraph->Draw("A*");
144  const string dimaName = translateWaveName(waveName);
145  cout << " drawing intensity for wave '" << waveName << "' -> '" << dimaName << "'" << endl;
146 
147  // draw Dima's plot
148  pad->cd();
149  TH1* hDima = dimaIntensityHists[dimaName];
150  if (!hDima) {
151  printWarn << "cannot find histogram with name '" << dimaName << "'. skipping." << endl;
152  return;
153  }
154  hDima->SetMarkerStyle(21);
155  hDima->SetMarkerSize(0.5);
156  hDima->SetLineColor(2);
157  hDima->SetMarkerColor(2);
158  hDima->Draw("SAME");
159 
160  TLegend *leg = new TLegend(0.82,0.97,0.97,0.82);
161  //leg->SetHeader("The Legend Title");
162  leg->AddEntry(rootgraph, "ROOTPWA","flpa");
163  leg->AddEntry(hDima,"COMPASSPWA","flpa");
164  leg->Draw();
165 
166  pad->Update();
167 
168  // readjust y-axis range
169  TMultiGraph* graph = static_cast<TMultiGraph*>(pad->GetPrimitive(waveName.c_str()));
170  if (graph) {
171  const double max = std::max(hDima->GetMaximum(), graph->GetHistogram()->GetMaximum());
172  graph->SetMaximum(1.1 * max);
173  pad->Modified();
174  pad->Update();
175  }
176 }
177 
178 void
179 compareIntensitiesWithDima(TTree* tree, // fitResult tree
180  TFile* dimaFile, const bool createPsFile)
181 {
182  printInfo << "comparing intensities with Dima's result." << endl;
183 
184  // get Dima's intensity histograms from file
185  map<string, TH1*> dimaIntensityHists; // intensity histograms w/o acceptance correction (see histgroups.txt)
186  {
187  map<string, TH1*> dimaHists = readDimaHists(dimaFile, TString("^h10\\d{2}$"));
188  // transform map with histogram titles as keys to map with wave names as keys
189  for (map<string, TH1*>::iterator i = dimaHists.begin(); i != dimaHists.end(); ++i) {
190  cout << " found intensity histogram '"<< i->first << "' " << flush;
191  TString waveName = i->first;
192  waveName = waveName.Strip(TString::kTrailing, ' ');
193  dimaIntensityHists[waveName.Data()] = i->second;
194  cout << "for wave '" << waveName << "'" << endl;
195  }
196  }
197 
198  // get ROOTPWA intensity histograms
199  vector<pair<string, TVirtualPad*> > wavePads = plotAllIntensities(tree, false, "./");
200 
201  // draw Dima's histograms on top of ROOTPWA's
202  cout << endl;
203  printInfo << "drawing Dima's histograms" << endl;
204  for (unsigned int i = 0; i < wavePads.size(); ++i) {
205  const string waveName = wavePads[i].first;
206  const string dimaName = translateWaveName(waveName);
207  cout << " drawing intensity for wave '" << waveName << "' -> '" << dimaName << "'" << endl;
208 
209  // draw Dima's plot
210  TLegend *leg = new TLegend(0.82,0.97,0.97,0.82);
211  TVirtualPad* pad = wavePads[i].second;
212  pad->cd();
213  TMultiGraph* list = (TMultiGraph*)pad->GetListOfPrimitives();
214  leg->AddEntry(&list[0], "ROOTPWA","flpa");
215  TH1* hDima = dimaIntensityHists[dimaName];
216  if (!hDima) {
217  printWarn << "cannot find histogram with name '" << dimaName << "'. skipping." << endl;
218  continue;
219  }
220  //hDima->Scale(3);
221  hDima->SetMarkerStyle(21);
222  hDima->SetMarkerSize(0.5);
223  hDima->SetLineColor(2);
224  hDima->SetMarkerColor(2);
225  hDima->Draw("SAME");
226 
227  leg->AddEntry(hDima,"COMPASSPWA","flpa");
228  leg->Draw();
229 
230  pad->Update();
231 
232  // readjust y-axis range
233  TMultiGraph* graph = static_cast<TMultiGraph*>(pad->GetPrimitive(waveName.c_str()));
234  if (graph) {
235  const double max = std::max(hDima->GetMaximum(), graph->GetHistogram()->GetMaximum());
236  graph->SetMaximum(1.1 * max);
237  pad->Modified();
238  pad->Update();
239  }
240  }
241  if (createPsFile) {
242  const string psFileName = "waveIntensities_compareDima.ps";
243  TCanvas dummyCanv("dummy", "dummy");
244  dummyCanv.Print((psFileName + "[").c_str());
245  vector<TCanvas*> canvases;
246  for (unsigned int i = 0; i < wavePads.size(); ++i) {
247  for (unsigned int j = 0; j < canvases.size(); j++) {
248  if(canvases[j] == wavePads[i].second->GetCanvas()) {
249  goto next;
250  }
251  }
252  canvases.push_back(wavePads[i].second->GetCanvas());
253  wavePads[i].second->GetCanvas()->Print(psFileName.c_str());
254  next: ;
255  }
256  dummyCanv.Print((psFileName + "]").c_str());
257  gSystem->Exec(("gv " + psFileName).c_str());
258  }
259 }
260 
261 
262 void
263 compareSpinTotalsWithDima(TTree* tree, // fitResult tree
264  TFile* dimaFile)
265 {
266  printInfo << "comparing spin totals with Dima's result." << endl;
267 
268  // get Dima's spin total histograms from file
269  map<string, TH1*> dimaSpinTotalHists;
270  {
271  map<string, TH1*> dimaHists = readDimaHists(dimaFile, TString("^h4\\d{3}$"));
272  // transform map with histogram titles as keys to map with wave names as keys
273  // this is still very simplistic and works only for part of the histograms
274  for (map<string, TH1*>::iterator i = dimaHists.begin(); i != dimaHists.end(); ++i) {
275  TString waveName = i->first;
276  bool waveMatch = true;
277  if (waveName == " Events ")
278  waveName = "";
279  else if ((waveName(1, 6) == "J^PC!=") && (waveName(18, 5) == "total")) {
280  waveName = waveName(8, 4);
281  waveName.ReplaceAll("^", "");
282  } else if ((waveName(1, 10) == "J^PC!M[c]=") && (waveName(22, 5) == "total")) {
283  waveName = waveName(12, 7);
284  waveName.ReplaceAll("^", "");
285  waveName.ReplaceAll("!", "");
286  } else
287  waveMatch = false;
288  if (waveMatch) {
289  cout << " found spin total histogram '"<< i->first << "' " << flush;
290  dimaSpinTotalHists[waveName.Data()] = i->second;
291  cout << "for wave '" << waveName << "'" << endl;
292  }
293  }
294  }
295 
296  // get ROOTPWA spinTotal histograms
297  vector<pair<string, TVirtualPad*> > wavePads = plotSpinTotals(tree, kBlack, 0, true, "");
298 
299  // draw Dima's histograms on top of ROOTPWA's
300  cout << endl;
301  printInfo<< "drawing Dima's histograms" << endl;
302  for (unsigned int i = 0; i < wavePads.size(); ++i) {
303  const string waveName = wavePads[i].first;
304  cout << " drawing spinTotal for wave '" << waveName << "'" << endl;
305 
306  // draw Dima's plot
307  TVirtualPad* pad = wavePads[i].second;
308  pad->cd();
309  TH1* hDima = dimaSpinTotalHists[waveName];
310  if (!hDima) {
311  printWarn << "cannot find histogram with name '" << waveName << "'. skipping." << endl;
312  continue;
313  }
314  hDima->Scale(3);
315  hDima->SetMarkerStyle(21);
316  hDima->SetMarkerSize(0.5);
317  hDima->SetLineColor(2);
318  hDima->SetMarkerColor(2);
319  hDima->Draw("SAME");
320  pad->Update();
321 
322  // readjust y-axis range
323  TString graphName = "total";
324  if (waveName != "") {
325  graphName = "g";
326  graphName.Append(waveName);
327  graphName.ReplaceAll("+", "p");
328  graphName.ReplaceAll("-", "m");
329  }
330  TGraphErrors* graph = static_cast<TGraphErrors*>(pad->GetPrimitive(graphName));
331  if (graph) {
332  const double max = std::max(hDima->GetMaximum(), graph->GetHistogram()->GetMaximum());
333  graph->SetMaximum(1.1 * max);
334  pad->Modified();
335  pad->Update();
336  }
337  }
338 }
339 
340 
341 void
342 comparePhasesWithDima(TTree* tree, // fitResult tree
343  TFile* dimaFile,
344  const vector<pair<string, string> >& wavenames_phases,
345  const string& branchName = "fitResult_v2")
346 {
347  printInfo << "comparing phases with Dima's result." << endl;
348 
349  // get Dima's phase histograms from file
350  map<string, TH1*> dimaPhaseHists; // relative phases (see histgroups.txt)
351  {
352  map<string, TH1*> dimaHists = readDimaHists(dimaFile, TString("^h4\\d{4}$"));
353  // transform map with histogram titles as keys to map with wave names as keys
354  for (map<string, TH1*>::iterator i = dimaHists.begin(); i != dimaHists.end(); ++i) {
355  cout << " found phase histogram '"<< i->first << "' " << flush;
356  TString hTitle = i->first;
357  hTitle = hTitle(4, hTitle.Length() - 5);
358  const int pos = hTitle.Index(" - ");
359  TString waveNames[2] = {hTitle(0, pos - 1), hTitle(pos + 3, hTitle.Length() - (pos + 3))};
360  waveNames[0] = waveNames[0].Strip(TString::kTrailing, ' ');
361  waveNames[1] = waveNames[1].Strip(TString::kTrailing, ' ');
362  dimaPhaseHists[(waveNames[0] + "|" + waveNames[1]).Data()] = i->second;
363  cout << "for waves '" << waveNames[0] << "', '" << waveNames[1] << "'" << endl;
364  }
365  }
366 
367  // get ROOTPWA phase histograms
368  fitResult* massBin = new fitResult();
369  tree->SetBranchAddress(branchName.c_str(), &massBin);
370  tree->GetEntry(0);
371 
372  for (unsigned int i = 0; i < wavenames_phases.size(); ++i) {
373  // get Dima's plot
374  const string histName = translateWaveName(wavenames_phases[i].first) + "|" + translateWaveName(wavenames_phases[i].second);
375  TH1* h = dimaPhaseHists[histName];
376  if (!h) {
377  printWarn << "cannot find histogram with name '" << histName << "'. skipping." << endl;
378  continue;
379  }
380 
381  // draw ROOTPWA phase histograms
382  new TCanvas();
383  TMultiGraph* rootgraph = plotPhase(tree, wavenames_phases[i].first, wavenames_phases[i].second);
384 
385  // draw Dima's plot on top of it
386  h->SetMarkerStyle(21);
387  h->SetMarkerSize(0.5);
388  h->SetLineColor(2);
389  h->SetMarkerColor(2);
390  h->Draw("SAME");
391 
392  TLegend *leg = new TLegend(0.82,0.97,0.97,0.82);
393  //leg->SetHeader("The Legend Title");
394  leg->AddEntry(rootgraph, "ROOTPWA","flpa");
395  leg->AddEntry(h,"COMPASSPWA","flpa");
396  leg->Draw();
397 
398  }
399 }
400 
401 
402 void
403 compareWithDima(TTree* tree, // fitResult tree
404  const string& dimaFileName = "/afs/e18.ph.tum.de/data/compass/bgrube/lambda/pwa/compassPWA/work/hfit.root",
405  const string& wavename = "",
406  const bool createPsFile =false)
407 {
408  if (!tree) {
409  printErr << "NULL pointer to tree. exiting." << endl;
410  return;
411  }
412 
413  printInfo << "opening Dima's file '" << dimaFileName << "'" << endl;
414  TFile* dimaFile = TFile::Open(dimaFileName.c_str(), "READ");
415  if (!dimaFile || dimaFile->IsZombie()) {
416  printErr << "cannot open file '" << dimaFileName << "'. exiting." << endl;
417  return;
418  }
419 
420  compareIntensityWithDima(tree, dimaFile, wavename);
421  //compareIntensitiesWithDima(tree, dimaFile, createPsFile);
422  //compareSpinTotalsWithDima(tree, dimaFile);
423  std::vector<std::pair<std::string, std::string> > phases;
424  phases.push_back(std::make_pair("1-1++0+rho770_01_pi0.amp", "1-2-+0+f21270_02_pi-.amp"));
425  phases.push_back(make_pair("1-2++1+rho770_21_pi0.amp", "1-1++0+rho770_01_pi0.amp"));
426  phases.push_back(make_pair("1-2++1+rho770_21_pi0.amp", "1-2-+0+f21270_02_pi-.amp"));
427  phases.push_back(make_pair("1-1-+1+rho770_11_pi0.amp", "1-2++1+rho770_21_pi0.amp"));
428  phases.push_back(make_pair("1-1-+1+rho770_11_pi0.amp", "1-1++0+rho770_01_pi0.amp"));
429 
430  comparePhasesWithDima(tree, dimaFile, phases);
431 
432 }
Generated by   doxygen 1.8.1.2