StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
compareRootFiles.cc
1 #include "compareRootFiles.hh"
2 
3 #include "TROOT.h"
4 #include "TStyle.h"
5 #include "TLatex.h"
6 #include "TLegend.h"
7 #include "TText.h"
8 #include "TAxis.h"
9 
10 #include <iostream>
11 #include <cstdlib>
12 
13 using std::cout;
14 using std::endl;
15 
16 compareRootFiles::compareRootFiles(string decayName, string newFileName, string oldFileName,
17  string description) {
18 
19  _decayName = decayName;
20  _newFileName = newFileName;
21  _oldFileName = oldFileName;
22  _description = description;
23 
24  _newFile = new TFile(newFileName.c_str(), "read");
25  _oldFile = new TFile(oldFileName.c_str(), "read");
26 
27  gROOT->SetStyle("Plain");
28  gStyle->SetOptStat(0);
29 
30  _theCanvas = new TCanvas("theCanvas", "", 900, 700);
31  _theCanvas->Clear();
32  _theCanvas->UseCurrentStyle();
33 
34  _emptyHist = new TH1D("", "", 10, 0.0, 0.0);
35 
36  _nGroupMax = 12;
37 
38 }
39 
40 compareRootFiles::~compareRootFiles() {
41 
42  _newFile->Close();
43  _oldFile->Close();
44 
45 }
46 
47 void compareRootFiles::getMtmPlots() {
48 
49  // Plot mtm plots for different particle groups
50  // Leptons+gamma, pions, kaons, charm, light/strange baryons, other
51 
52  vector<TH1D*> newHistVect, oldHistVect;
53 
54  vector<int> leptonInts;
55  leptonInts.push_back(0); leptonInts.push_back(1);
56  TH1D* newLeptMtmHist = getMtmHist(_newFile, "newLeptMtmHist", leptonInts);
57  TH1D* oldLeptMtmHist = getMtmHist(_oldFile, "oldLeptMtmHist", leptonInts);
58  newHistVect.push_back(newLeptMtmHist);
59  oldHistVect.push_back(oldLeptMtmHist);
60 
61  vector<int> pionInts;
62  pionInts.push_back(2); pionInts.push_back(3);
63  TH1D* newPionMtmHist = getMtmHist(_newFile, "newPionMtmHist", pionInts);
64  TH1D* oldPionMtmHist = getMtmHist(_oldFile, "oldPionMtmHist", pionInts);
65  newHistVect.push_back(newPionMtmHist);
66  oldHistVect.push_back(oldPionMtmHist);
67 
68  vector<int> kaonInts;
69  kaonInts.push_back(4); kaonInts.push_back(5);
70  TH1D* newKaonMtmHist = getMtmHist(_newFile, "newKaonMtmHist", kaonInts);
71  TH1D* oldKaonMtmHist = getMtmHist(_oldFile, "oldKaonMtmHist", kaonInts);
72  newHistVect.push_back(newKaonMtmHist);
73  oldHistVect.push_back(oldKaonMtmHist);
74 
75  vector<int> charmInts;
76  charmInts.push_back(6); charmInts.push_back(7);
77  TH1D* newCharmMtmHist = getMtmHist(_newFile, "newCharmMtmHist", charmInts);
78  TH1D* oldCharmMtmHist = getMtmHist(_oldFile, "oldCharmMtmHist", charmInts);
79  newHistVect.push_back(newCharmMtmHist);
80  oldHistVect.push_back(oldCharmMtmHist);
81 
82  vector<int> baryonInts;
83  baryonInts.push_back(8); baryonInts.push_back(9);
84  TH1D* newBaryonMtmHist = getMtmHist(_newFile, "newBaryonMtmHist", baryonInts);
85  TH1D* oldBaryonMtmHist = getMtmHist(_oldFile, "oldBaryonMtmHist", baryonInts);
86  newHistVect.push_back(newBaryonMtmHist);
87  oldHistVect.push_back(oldBaryonMtmHist);
88 
89  vector<int> otherInts;
90  otherInts.push_back(10);
91  TH1D* newOtherMtmHist = getMtmHist(_newFile, "newOtherMtmHist", otherInts);
92  TH1D* oldOtherMtmHist = getMtmHist(_oldFile, "oldOtherMtmHist", otherInts);
93  newHistVect.push_back(newOtherMtmHist);
94  oldHistVect.push_back(oldOtherMtmHist);
95 
96  _theCanvas->Clear();
97  _theCanvas->UseCurrentStyle();
98  _theCanvas->cd(1);
99  gPad->SetLogy();
100 
101  // Different colours for the particle groups
102  int colours[6] = {1, 2, 4, 6, 8, 28};
103  string histLabels[6] = {"Leptons,#gamma", "Pions", "Kaons", "D+,D-,D0",
104  "Light/strange baryons", "Other"};
105 
106  TLegend* theLegend = new TLegend(0.7, 0.7, 0.9, 0.9);
107  theLegend->SetFillColor(kWhite);
108 
109  double newHistMax(0.0), oldHistMax(0.0);
110  int newHistMaxInt(0), oldHistMaxInt(0);
111 
112  double totalNewFreq(0.0), totalOldFreq(0.0);
113 
114  int iHist(0);
115  int nHistVect = newHistVect.size();
116  for (iHist = 0; iHist < nHistVect; iHist++) {
117 
118  TH1D* newHist = newHistVect[iHist];
119  string labelName = histLabels[iHist];
120  theLegend->AddEntry(newHist, labelName.c_str(), "lpf");
121 
122  int colour = colours[iHist];
123  newHist->SetLineColor(colour);
124 
125  double theNewHistMax = newHist->GetMaximum();
126  if (theNewHistMax > newHistMax) {
127  newHistMax = theNewHistMax;
128  newHistMaxInt = iHist;
129  }
130 
131  totalNewFreq += newHist->Integral();
132 
133  TH1D* oldHist = oldHistVect[iHist];
134  oldHist->SetLineStyle(kDashed);
135  oldHist->SetLineColor(colour);
136 
137  double theOldHistMax = oldHist->GetMaximum();
138  if (theOldHistMax > oldHistMax) {
139  oldHistMax = theOldHistMax;
140  oldHistMaxInt = iHist;
141  }
142 
143  totalOldFreq += oldHist->Integral();
144 
145  }
146 
147  double newScale = 1.0/totalNewFreq;
148  double oldScale = 1.0/totalOldFreq;
149 
150  if (newHistMax > oldHistMax) {
151 
152  TH1D* newFirstHist = newHistVect[newHistMaxInt];
153 
154  TAxis* xAxis = newFirstHist->GetXaxis();
155  double dx = xAxis->GetBinWidth(1)*1000.0;
156  char dxChar[100];
157  sprintf(dxChar, "Normalised frequency/%.0f (MeV/c)", dx);
158 
159  newFirstHist->SetXTitle("p (GeV/c)");
160  newFirstHist->SetYTitle(dxChar);
161  newFirstHist->SetTitleOffset(1.25, "Y");
162  newFirstHist->Scale(newScale);
163  newFirstHist->Draw();
164 
165  for (iHist = 0; iHist < nHistVect; iHist++) {
166 
167  if (iHist != newHistMaxInt) {
168  TH1D* otherNewHist = newHistVect[iHist];
169  otherNewHist->Scale(newScale);
170  otherNewHist->Draw("same");
171 
172  }
173 
174  TH1D *otherOldHist = oldHistVect[iHist];
175  otherOldHist->Scale(oldScale);
176  otherOldHist->Draw("same");
177 
178  }
179 
180  } else {
181 
182  TH1D* oldFirstHist = oldHistVect[oldHistMaxInt];
183 
184  TAxis* xAxis = oldFirstHist->GetXaxis();
185  double dx = xAxis->GetBinWidth(1)*1000.0;
186  char dxChar[100];
187  sprintf(dxChar, "Normalised frequency/%.0f (MeV/c)", dx);
188 
189  oldFirstHist->SetXTitle("p (GeV/c)");
190  oldFirstHist->SetYTitle(dxChar);
191  oldFirstHist->SetTitleOffset(1.25, "Y");
192  oldFirstHist->Scale(oldScale);
193  oldFirstHist->Draw();
194 
195  for (iHist = 0; iHist < nHistVect; iHist++) {
196 
197  if (iHist != oldHistMaxInt) {
198  TH1D* otherOldHist = oldHistVect[iHist];
199  otherOldHist->Scale(oldScale);
200  otherOldHist->Draw("same");
201  }
202 
203  TH1D* otherNewHist = newHistVect[iHist];
204  otherNewHist->Scale(newScale);
205  otherNewHist->Draw("same");
206 
207  }
208 
209  }
210 
211  theLegend->Draw();
212 
213  TText* text = new TText();
214  text->SetTextSize(0.03);
215  text->DrawTextNDC(0.1, 0.95, "New version = solid lines");
216  text->DrawTextNDC(0.1, 0.915, "Old version = dotted lines");
217  TLatex* latexString = new TLatex();
218  latexString->SetTextSize(0.03);
219  latexString->SetNDC();
220  latexString->DrawLatex(0.5, 0.92, _description.c_str());
221 
222  string gifFileName("gifFiles/");
223  gifFileName.append(_decayName); gifFileName.append("_mtmHist.gif");
224  _theCanvas->Print(gifFileName.c_str());
225 
226  gPad->SetLogy(0);
227 
228 }
229 
230 void compareRootFiles::getPartIdPlots() {
231 
232  TH1D* newPartIdHist = getPartIdHist(_newFile, "newPartIdHist");
233  TH1D* oldPartIdHist = getPartIdHist(_oldFile, "oldPartIdHist");
234 
235  _theCanvas->Clear();
236  _theCanvas->UseCurrentStyle();
237  _theCanvas->cd(1);
238  gPad->SetLogy();
239 
240  double newPartIdMax = newPartIdHist->GetMaximum();
241  double oldPartIdMax = oldPartIdHist->GetMaximum();
242 
243  oldPartIdHist->SetLineStyle(kDashed);
244 
245  if (newPartIdMax > oldPartIdMax) {
246 
247  newPartIdHist->SetXTitle("");
248  newPartIdHist->SetYTitle("Normalised frequency");
249  newPartIdHist->SetTitleOffset(1.25, "Y");
250 
251  newPartIdHist->Draw();
252  oldPartIdHist->Draw("same");
253 
254  } else {
255 
256  oldPartIdHist->SetXTitle("");
257  oldPartIdHist->SetYTitle("Normalised frequency");
258  oldPartIdHist->SetTitleOffset(1.25, "Y");
259 
260  oldPartIdHist->Draw();
261  newPartIdHist->Draw("same");
262 
263  }
264 
265  TLatex* latex = new TLatex();
266  latex->SetTextAngle(90.0);
267 
268  latex->SetTextSize(0.035);
269  latex->DrawLatex(0.5, 0, "Leptons");
270 
271  latex->SetTextSize(0.045);
272  latex->DrawLatex(1.5, 0, "#gamma");
273  latex->DrawLatex(2.5, 0, "#pi^{#pm}");
274  latex->DrawLatex(3.5, 0, "#pi^{0}");
275  latex->DrawLatex(4.5, 0, "K^{#pm}");
276  latex->DrawLatex(5.5, 0, "K^{0}");
277  latex->DrawLatex(6.5, 0, "D^{#pm}");
278  latex->DrawLatex(7.5, 0, "D^{0}");
279 
280  latex->SetTextSize(0.035);
281  latex->DrawLatex(8.4, 0, "Light"); latex->DrawLatex(8.75, 0, "baryons");
282  latex->DrawLatex(9.4, 0, "Strange"); latex->DrawLatex(9.75, 0, "baryons");
283 
284  latex->DrawLatex(10.5, 0, "Other");
285 
286  TText* text = new TText();
287  text->SetTextSize(0.03);
288  text->DrawTextNDC(0.1, 0.95, "New version = solid lines");
289  text->DrawTextNDC(0.1, 0.915, "Old version = dotted lines");
290  TLatex* latexString = new TLatex();
291  latexString->SetTextSize(0.03);
292  latexString->SetNDC();
293  latexString->DrawLatex(0.5, 0.92, _description.c_str());
294 
295  string gifFileName("gifFiles/");
296  gifFileName.append(_decayName); gifFileName.append("_partIdHist.gif");
297  _theCanvas->Print(gifFileName.c_str());
298 
299  gPad->SetLogy(0);
300 
301 }
302 
303 void compareRootFiles::getAllIdPlots() {
304 
305  TH1D* newAllIdHist = getAllIdHist(_newFile, "newAllIdHist");
306  TH1D* oldAllIdHist = getAllIdHist(_oldFile, "oldAllIdHist");
307 
308  _theCanvas->Clear();
309  _theCanvas->UseCurrentStyle();
310  _theCanvas->cd(1);
311  gPad->SetLogy();
312 
313  double newAllIdMax = newAllIdHist->GetMaximum();
314  double oldAllIdMax = oldAllIdHist->GetMaximum();
315 
316  oldAllIdHist->SetLineStyle(kDashed);
317 
318  if (newAllIdMax > oldAllIdMax) {
319 
320  newAllIdHist->SetXTitle("PDG Id");
321  newAllIdHist->SetYTitle("Normalised frequency");
322  newAllIdHist->SetTitleOffset(1.25, "Y");
323 
324  newAllIdHist->Draw();
325  oldAllIdHist->Draw("same");
326 
327  } else {
328 
329  oldAllIdHist->SetXTitle("PDG Id");
330  oldAllIdHist->SetYTitle("Normalised frequency");
331  oldAllIdHist->SetTitleOffset(1.25, "Y");
332 
333  oldAllIdHist->Draw();
334  newAllIdHist->Draw("same");
335 
336  }
337 
338  TText* text = new TText();
339  text->SetTextSize(0.03);
340  text->DrawTextNDC(0.1, 0.95, "New version = solid lines");
341  text->DrawTextNDC(0.1, 0.915, "Old version = dotted lines");
342  TLatex* latexString = new TLatex();
343  latexString->SetTextSize(0.03);
344  latexString->SetNDC();
345  latexString->DrawLatex(0.5, 0.92, _description.c_str());
346 
347  string gifFileName("gifFiles/");
348  gifFileName.append(_decayName); gifFileName.append("_allIdHist.gif");
349  _theCanvas->Print(gifFileName.c_str());
350 
351  gPad->SetLogy(0);
352 
353 }
354 
355 void compareRootFiles::getNDaugPlots() {
356 
357  TH1D* newDaugHist = getDaugHist(_newFile, "newDaugHist");
358  TH1D* oldDaugHist = getDaugHist(_oldFile, "oldDaugHist");
359 
360  _theCanvas->Clear();
361  _theCanvas->UseCurrentStyle();
362  _theCanvas->cd(1);
363 
364  double newDaugMax = newDaugHist->GetMaximum();
365  double oldDaugMax = oldDaugHist->GetMaximum();
366 
367  oldDaugHist->SetLineStyle(kDashed);
368 
369  if (newDaugMax > oldDaugMax) {
370 
371  newDaugHist->SetXTitle("nDaughters");
372  newDaugHist->SetYTitle("Normalised frequency");
373  newDaugHist->SetTitleOffset(1.25, "Y");
374 
375  newDaugHist->Draw();
376  oldDaugHist->Draw("same");
377 
378  } else {
379 
380  oldDaugHist->SetXTitle("nDaughters");
381  oldDaugHist->SetYTitle("Normalised frequency");
382  oldDaugHist->SetTitleOffset(1.25, "Y");
383 
384  oldDaugHist->Draw();
385  newDaugHist->Draw("same");
386 
387  }
388 
389  TText* text = new TText();
390  text->SetTextSize(0.03);
391  text->DrawTextNDC(0.1, 0.95, "New version = solid lines");
392  text->DrawTextNDC(0.1, 0.915, "Old version = dotted lines");
393  TLatex* latexString = new TLatex();
394  latexString->SetTextSize(0.03);
395  latexString->SetNDC();
396  latexString->DrawLatex(0.5, 0.92, _description.c_str());
397 
398  string gifFileName("gifFiles/");
399  gifFileName.append(_decayName); gifFileName.append("_daugHist.gif");
400  _theCanvas->Print(gifFileName.c_str());
401 
402 }
403 
404 TH1D* compareRootFiles::getMtmHist(TFile* theFile, string histName, vector<int> groupInts) {
405 
406  if (theFile == 0) {
407  // Return empty histogram
408  return _emptyHist;
409  }
410 
411  TTree* theTree = dynamic_cast<TTree*>(theFile->Get("Data"));
412  if (theTree == 0) {
413  // Return empty histogram
414  return _emptyHist;
415  }
416 
417  int id;
418  double p;
419  theTree->SetBranchAddress("id", &id);
420  theTree->SetBranchAddress("p", &p);
421 
422  double pMax = theTree->GetMaximum("p");
423  TH1D* mtmHist = new TH1D(histName.c_str(), "", 100, 0.0, pMax*1.1);
424 
425  int i;
426  int nEntries = theTree->GetEntries();
427  for (i = 0; i < nEntries; i++) {
428 
429  theTree->GetEntry(i);
430  int testInt = this->getPartGroup(id);
431 
432  // See if the particle id matches any in the group integer vector
433  vector<int>::iterator iter;
434  for (iter = groupInts.begin(); iter != groupInts.end(); ++iter) {
435  int groupInt = *iter;
436 
437  if (groupInt == testInt) {
438 
439  // We have the right particle group id code
440  mtmHist->Fill(p);
441  }
442  }
443 
444  }
445 
446  return mtmHist;
447 
448 }
449 
450 TH1D* compareRootFiles::getPartIdHist(TFile* theFile, string histName) {
451 
452  // Group pi,K,D etc.. into groups and plot multiplicities
453  // with bins along the x axis representing each group
454  // leptons gamma pi+- pi0 K+- K0 D+- D0 light_baryons strange_baryons other
455 
456  if (theFile == 0) {
457  // Return empty histogram
458  return _emptyHist;
459  }
460 
461  TTree* theTree = dynamic_cast<TTree*>(theFile->Get("Data"));
462  if (theTree == 0) {
463  // Return empty histogram
464  return _emptyHist;
465  }
466 
467  int id;
468  theTree->SetBranchAddress("id", &id);
469 
470  int nGroupMax = _nGroupMax;
471  vector<double> partNumbers(nGroupMax);
472  int iP;
473  for (iP = 0; iP < nGroupMax; iP++) {
474  partNumbers[iP] = 0.0;
475  }
476 
477  int i;
478  int nEntries = theTree->GetEntries();
479  for (i = 0; i < nEntries; i++) {
480 
481  theTree->GetEntry(i);
482 
483  int groupInt = this->getPartGroup(id);
484  if (groupInt >= 0 && groupInt < nGroupMax) {
485  partNumbers[groupInt] += 1.0;
486  }
487 
488  }
489 
490  TH1D* idHist = new TH1D(histName.c_str(), "", nGroupMax, 0, nGroupMax);
491  idHist->SetMinimum(1);
492 
493  for (iP = 0; iP < nGroupMax; iP++) {
494 
495  double total = partNumbers[iP];
496  idHist->Fill(iP, total);
497 
498  }
499 
500  idHist->Scale(1.0/(nEntries*1.0));
501  idHist->SetMaximum(1.0);
502 
503  return idHist;
504 
505 }
506 
507 int compareRootFiles::getPartGroup(int PDGId) {
508 
509  int group(-1);
510 
511  int absPDGId = std::abs(PDGId);
512 
513  if (absPDGId >= 11 && absPDGId <= 16) {
514  group = 0; // leptons
515  } else if (absPDGId == 22) {
516  group = 1; // photon
517  } else if (absPDGId == 211) {
518  group = 2; // pi+-
519  } else if (absPDGId == 111) {
520  group = 3; // pi0
521  } else if (absPDGId == 321) {
522  group = 4; // K+-
523  } else if (absPDGId == 311 || absPDGId == 130 || absPDGId == 310) {
524  group = 5; // K0
525  } else if (absPDGId == 411) {
526  group = 6; // D+-
527  } else if (absPDGId == 421) {
528  group = 7; // D0
529  } else if (absPDGId == 2212 || absPDGId == 2112 || absPDGId == 2224 ||
530  absPDGId == 2214 || absPDGId == 2114 || absPDGId == 1114) {
531  group = 8; // light baryons
532  } else if (absPDGId >= 3112 && absPDGId <= 3334) {
533  group = 9; // strange baryons
534  } else if (absPDGId != 0) {
535  group = 10; // other particles
536  }
537 
538  return group;
539 
540 }
541 
542 TH1D* compareRootFiles::getAllIdHist(TFile* theFile, string histName) {
543 
544  if (theFile == 0) {
545  // Return empty histogram
546  return _emptyHist;
547  }
548 
549  TTree* theTree = dynamic_cast<TTree*>(theFile->Get("Data"));
550  if (theTree == 0) {
551  // Return empty histogram
552  return _emptyHist;
553  }
554 
555  int id;
556  theTree->SetBranchAddress("id", &id);
557 
558  // Just limit the PDG id integer codes up to 4000
559  int idMax(4000);
560  TH1D* idHist = new TH1D(histName.c_str(), "", 100, -idMax, idMax);
561 
562  int nEntries = theTree->GetEntries();
563  cout<<"Number of entries = "<<nEntries<<endl;
564  int i;
565  for (i = 0; i < nEntries; i++) {
566 
567  theTree->GetEntry(i);
568  idHist->Fill(id);
569 
570  }
571 
572  double nIdScale = idHist->Integral();
573  idHist->Scale(1.0/nIdScale);
574  idHist->SetMaximum(1.0);
575 
576  return idHist;
577 
578 }
579 
580 TH1D* compareRootFiles::getDaugHist(TFile* theFile, string histName) {
581 
582  if (theFile == 0) {
583  // Return empty histogram
584  return _emptyHist;
585  }
586 
587  TTree* nDaugTree = dynamic_cast<TTree*>(theFile->Get("nDaugTree"));
588  if (nDaugTree == 0) {
589  // Return empty histogram
590  return _emptyHist;
591  }
592 
593  int nDaug;
594  nDaugTree->SetBranchAddress("nDaug", &nDaug);
595 
596  int nDaugMax = (int) nDaugTree->GetMaximum("nDaug");
597  int nDaugLimit = nDaugMax + 2;
598 
599  TH1D* nDaugHist = new TH1D(histName.c_str(), "", nDaugLimit, 0, nDaugLimit);
600 
601  int nEntries = nDaugTree->GetEntries();
602  cout<<"Number of entries = "<<nEntries<<endl;
603  int i;
604  for (i = 0; i < nEntries; i++) {
605 
606  nDaugTree->GetEntry(i);
607 
608  if (nDaug > 0) {
609  nDaugHist->Fill(nDaug);
610  }
611 
612  }
613 
614  double nDaugScale = nDaugHist->Integral();
615  nDaugHist->Scale(1.0/nDaugScale);
616 
617  cout<<"nDaugScale = "<<nDaugScale<<endl;
618 
619  return nDaugHist;
620 }
621 
622 
623 int main(int argc, char** argv) {
624 
625  string decayName("UpsilonDecay1");
626  if (argc > 1) {decayName = argv[1];}
627 
628  string newRootFile("rootFiles/newUpsilonDecay1.root");
629  if (argc > 2) {newRootFile = argv[2];}
630 
631  string oldRootFile("rootFiles/oldUpsilonDecay1.root");
632  if (argc > 3) {oldRootFile = argv[3];}
633 
634  string description("Example description");
635  if (argc > 4) {description = argv[4];}
636 
637  compareRootFiles myCompareRootFiles(decayName, newRootFile, oldRootFile, description);
638 
639  myCompareRootFiles.getNDaugPlots();
640  myCompareRootFiles.getAllIdPlots();
641  myCompareRootFiles.getPartIdPlots();
642  myCompareRootFiles.getMtmPlots();
643 
644 }
Generated by   doxygen 1.8.5