File indexing completed on 2024-04-06 12:33:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "TH1F.h"
0011 #include "TProfile.h"
0012 #include "TCanvas.h"
0013 #include "TMath.h"
0014 #include "TPad.h"
0015 #include "TPaveLabel.h"
0016 #include <vector>
0017
0018 void NormalizeHistogramsToFirst(TH1* h1, TH1* h2);
0019 void NormalizeHistogramsTo1(TH1* h1, TH1* h2);
0020 TH1* PlotRatiosHistograms(TH1* h1, TH1* h2);
0021
0022 int ratioCounter = 0;
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036 void SetGlobalStyle() {
0037 gROOT->SetStyle("Plain");
0038 gStyle->SetPadGridX(kTRUE);
0039 gStyle->SetPadGridY(kTRUE);
0040 gStyle->SetPadRightMargin(0.07);
0041 gStyle->SetPadLeftMargin(0.13);
0042
0043
0044
0045
0046
0047 gStyle->SetOptStat(0);
0048 }
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059 void SetHistogramStyle(TH1* h, Style_t mstyle, Color_t color, Size_t msize = 0.7, Width_t lwidth = 2,
0060 Float_t tsize = 0.05, Float_t toffset = 1.2) {
0061 if (!h)
0062 return;
0063 h->SetMarkerStyle(mstyle);
0064 h->SetMarkerColor(color);
0065 h->SetMarkerSize(msize);
0066 h->SetLineColor(color);
0067 h->SetLineWidth(lwidth);
0068 h->GetYaxis()->SetTitleSize(tsize);
0069 h->GetYaxis()->SetTitleOffset(toffset);
0070 h->GetXaxis()->SetLabelFont(63);
0071 h->GetXaxis()->SetLabelSize(14);
0072 h->GetYaxis()->SetLabelFont(63);
0073 h->GetYaxis()->SetLabelSize(14);
0074 }
0075
0076
0077
0078
0079
0080
0081
0082 TList* GetListOfDirectories(TDirectory* dir, const char* match = 0) {
0083 TIter nextkey(dir->GetListOfKeys());
0084 TList* sl = new TList();
0085 TKey* key = 0;
0086 TKey* oldkey = 0;
0087
0088 while (( key = (TKey*)nextkey() )) {
0089 TObject *obj = key->ReadObj();
0090 if ( obj->IsA()->InheritsFrom( "TDirectory" ) ) {
0091 TString theName = obj->GetName();
0092 if (!match) {
0093 cout << " -> " << theName << endl;
0094 sl->Add(obj);
0095 }
0096 else if (theName.Contains(match)) {
0097 cout << " -> " << theName << endl;
0098 sl->Add(obj);
0099 }
0100 }
0101 }
0102
0103 return sl;
0104 }
0105
0106
0107
0108
0109
0110
0111
0112 TList* getListOfBranches(const char* dataType, TFile* file, const char* branchContent) {
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140 if (TString(dataType) == "RECO") {
0141 if(! file->cd("DQMData/Run 1/Muons/Run summary")) {
0142 cout << "ERROR: Muon Histos for " << dataType << " not found" << endl;
0143 return 0;
0144 }
0145 }
0146 else {
0147 cout << "ERROR: Data type " << dataType << " not allowed: only RECO is considered" << endl;
0148 return 0;
0149 }
0150
0151 TDirectory * dir=gDirectory;
0152 TList* sl = GetListOfDirectories(dir, branchContent);
0153
0154 if (sl->GetSize() == 0) {
0155 cout << "ERROR: No DQM muon reco histos found in NEW file " << endl;
0156 delete sl;
0157 return 0;
0158 }
0159
0160 return sl;
0161 }
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173 double KolmogorovTest(TH1 *h1, TH1 *h2){
0174
0175 double mya_array[1300], myb_array[1300];
0176 vector<double> mya;
0177 vector<double> myb;
0178
0179
0180 for (int i=0; i<h1->GetNbinsX(); i++){
0181 mya.push_back(h1->GetBinContent(i+1));
0182 myb.push_back(h2->GetBinContent(i+1));
0183 }
0184
0185 sort(mya.begin(),mya.end());
0186 sort(myb.begin(),myb.end());
0187 copy(mya.begin(),mya.end(),mya_array);
0188 copy(myb.begin(),myb.end(),myb_array);
0189
0190 const int nbinsa = h1->GetNbinsX();
0191 const int nbinsb = h2->GetNbinsX();
0192
0193 double kstest = TMath::KolmogorovTest(nbinsa, mya_array,
0194 nbinsb, myb_array,
0195 "UOX");
0196 #ifdef DEBUG
0197 cout << "DEBUG: + KS value = " << kstest << endl;
0198 #endif
0199
0200
0201
0202 TString legend = Form("KS=%4.2f", kstest);
0203
0204
0205
0206 TPaveLabel* pl = new TPaveLabel(0.79,0.91,0.93,0.96, legend.Data(), "NDC");
0207
0208
0209
0210 pl->SetLineColor(41);
0211 pl->SetLineWidth(1);
0212 pl->SetLineStyle(1);
0213 pl->SetFillColor(41);
0214 pl->SetBorderSize(3);
0215
0216 if (kstest < 0.7)
0217 pl->SetTextColor(kRed);
0218
0219 pl->Draw();
0220
0221 return kstest;
0222 }
0223
0224
0225
0226
0227
0228
0229
0230 void setStats(TH1* s, TH1* r, double startingY, double startingX = .1, bool fit = false){
0231 if (startingY<0){
0232 s->SetStats(0);
0233 r->SetStats(0);
0234 }
0235 else {
0236
0237
0238 if (fit){
0239 s->Fit("gaus");
0240 TF1* f1 = (TF1*) s->GetListOfFunctions()->FindObject("gaus");
0241 if (f1) {
0242 f1->SetLineColor(2);
0243 f1->SetLineWidth(1);
0244 }
0245 }
0246 s->Draw();
0247 gPad->Update();
0248 TPaveStats* st1 = (TPaveStats*) s->GetListOfFunctions()->FindObject("stats");
0249 if (st1) {
0250 if (fit) {st1->SetOptFit(0010); st1->SetOptStat(1001);}
0251 st1->SetX1NDC(startingX);
0252 st1->SetX2NDC(startingX+0.30);
0253 st1->SetY1NDC(startingY+0.20);
0254 st1->SetY2NDC(startingY+0.35);
0255 st1->SetTextColor(2);
0256 }
0257 else s->SetStats(0);
0258 if (fit) {
0259 r->Fit("gaus");
0260 TF1* f2 = (TF1*) r->GetListOfFunctions()->FindObject("gaus");
0261 if (f2) {
0262 f2->SetLineColor(4);
0263 f2->SetLineWidth(1);
0264 }
0265 }
0266 r->Draw();
0267 gPad->Update();
0268 TPaveStats* st2 = (TPaveStats*) r->GetListOfFunctions()->FindObject("stats");
0269 if (st2) {
0270 if (fit) {st2->SetOptFit(0010); st2->SetOptStat(1001);}
0271 st2->SetX1NDC(startingX);
0272 st2->SetX2NDC(startingX+0.30);
0273 st2->SetY1NDC(startingY);
0274 st2->SetY2NDC(startingY+0.15);
0275 st2->SetTextColor(4);
0276 }
0277 else r->SetStats(0);
0278 }
0279 }
0280
0281
0282
0283
0284
0285
0286
0287 void PlotNHistograms(const TString& pdfFile,
0288 TDirectory* rdir, TDirectory* sdir,
0289 const TString& rcollname, const TString& scollname,
0290 const char* canvasName, const char* canvasTitle,
0291 const TString& refLabel, const TString& newLabel,
0292 unsigned int nhistos, const char** hnames, const char** htitles = 0,
0293 bool* logy = 0, bool* doKolmo = 0, Double_t* norm = 0, bool *resol = 0,
0294 Double_t* minx = 0, Double_t* maxx = 0,
0295 Double_t* miny = 0, Double_t* maxy = 0) {
0296 #ifdef DEBUG
0297 cout << " + Plotting histograms for " << canvasTitle << endl;
0298 cerr << " + Plotting histograms for " << canvasTitle << endl;
0299 #endif
0300 TH1* rh = 0;
0301 TH1* sh = 0;
0302
0303
0304
0305
0306
0307 #ifdef DEBUG
0308 cout << "DEBUG: Creating canvas..." << endl;
0309 #endif
0310 TCanvas* canvas = 0;
0311 if (nhistos >4)
0312 canvas = new TCanvas(canvasName, canvasTitle, 1000, 1400);
0313 else
0314 canvas = new TCanvas(canvasName, canvasTitle, 1000, 1050);
0315 canvas->Divide(2,(nhistos+1)/2);
0316
0317 #ifdef DEBUG
0318 cout << "DEBUG: Looping over histograms" << endl;
0319 #endif
0320 for (unsigned int i = 0; i < nhistos; i++) {
0321 #ifdef DEBUG
0322 cout << "DEBUG: [" << i << "] histogram name: " << flush << hnames[i] << endl;
0323 cout << "DEBUG: rdir: " << rdir << endl;
0324 cout << "DEBUG: sdir: " << sdir << endl;
0325 cout << "DEBUG: rcollname: " << rcollname << endl;
0326 cout << "DEBUG: scollname: " << scollname << endl;
0327 cout << "DEBUG: rname: " << (rcollname + "/" + hnames[i]) << endl;
0328 cout << "DEBUG: sname: " << (scollname + "/" + hnames[i]) << endl;
0329 #endif
0330
0331 TString hnamesi = hnames[i];
0332
0333
0334 if (hnamesi == "") continue;
0335
0336
0337
0338 TString hnames_tmp = hnamesi;
0339 rdir->cd(rcollname);
0340 TIter next(gDirectory->GetListOfKeys());
0341 TKey *key;
0342 while ((key = (TKey*)next())) {
0343
0344 TObject *obj = key->ReadObj();
0345 TString temp = obj->GetName();
0346 if ( (temp.Contains("nhits_vs_eta_pfx")) &&
0347 hnamesi.Contains("hits_eta")
0348 ) {
0349 hnames_tmp = temp;
0350 }
0351 if ( (temp.Contains("chi2_vs_eta_pfx")) &&
0352 hnamesi.Contains("chi2mean")
0353 ) {
0354 hnames_tmp = temp;
0355 }
0356 }
0357
0358
0359
0360 #ifdef DEBUG
0361 cout << "DEBUG: Getting object reference samples " << (rcollname + "/" + hnames_tmp) << " from " << rdir << endl;
0362 #endif
0363 rdir->GetObject(rcollname + "/" + hnames_tmp, rh);
0364 if (! rh) {
0365 cout << "WARNING: Could not find a reference histogram or profile named " << hnames_tmp
0366 << " in " << rdir->GetName() << endl;
0367 cout << " Skipping" << endl;
0368 continue;
0369 }
0370
0371
0372 if (TString(rh->IsA()->GetName()) == "TH2F") {
0373 #ifdef DEBUG
0374 cout << "DEBUG: It is a TH2F object... project in Y!" << endl;
0375 #endif
0376 TH1* proj = ((TH2F*) rh)->ProjectionY();
0377 rh = proj;
0378 }
0379
0380
0381
0382
0383 hnames_tmp=hnamesi;
0384
0385 sdir->cd(scollname);
0386 TIter next2(gDirectory->GetListOfKeys());
0387 TKey *key2;
0388 while ((key2 = (TKey*)next2())) {
0389
0390 TObject *obj = key2->ReadObj();
0391 TString temp = obj->GetName();
0392 if ( (temp.Contains("nhits_vs_eta_pfx")) &&
0393 hnamesi.Contains("hits_eta")
0394 ) {
0395 hnames_tmp = temp;
0396 }
0397 if ( (temp.Contains("chi2_vs_eta_pfx")) &&
0398 hnamesi.Contains("chi2mean")
0399 ) {
0400 hnames_tmp = temp;
0401 }
0402 }
0403
0404 #ifdef DEBUG
0405 cout << "DEBUG: Getting object for selected sample " << (scollname + "/" + hnames_tmp)
0406 << " from " << sdir << endl;
0407 #endif
0408 sdir->GetObject(scollname + "/" + hnames_tmp, sh);
0409 if (! sh) {
0410 cout << "WARNING: Could not find a signal histogram or profile named " << hnames_tmp
0411 << " in " << rdir->GetName() << endl;
0412 cout << " Skipping" << endl;
0413 continue;
0414 }
0415
0416
0417 if (TString(sh->IsA()->GetName()) == "TH2F") {
0418 #ifdef DEBUG
0419 cout << "DEBUG: " << hnames_tmp << " is a TH2F object... project in Y!" << endl;
0420 #endif
0421 TH1* proj = ((TH2F*) sh)->ProjectionY();
0422 sh = proj;
0423 }
0424
0425 if(TString(sh->GetName()).Contains(" vs "))norm[i]= -999.;
0426
0427
0428
0429 if (norm[i] == -1.)
0430 NormalizeHistogramsTo1(rh, sh);
0431 else if (norm[i] == 0.)
0432 NormalizeHistogramsToFirst(rh,sh);
0433 else if (norm[i] == -999.){
0434 cout << "DEBUG: Normalizing histograms to nothing" << "..." << endl;
0435 }
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446 #ifdef DEBUG
0447 cout << "DEBUG: Setting styles..." << endl;
0448 #endif
0449 SetHistogramStyle(rh, 21, 4);
0450 SetHistogramStyle(sh, 20, 2);
0451
0452 if (htitles) {
0453 rh->SetTitle(htitles[i]);
0454 sh->SetTitle(htitles[i]);
0455 rh->GetYaxis()->SetTitle(htitles[i]);
0456 sh->GetYaxis()->SetTitle(htitles[i]);
0457 }
0458
0459 if (minx) {
0460
0461
0462
0463 rh->GetXaxis()->SetRangeUser(minx[i],maxx[i]);
0464 sh->GetXaxis()->SetRangeUser(minx[i],maxx[i]);
0465
0466
0467 }
0468
0469 if (miny) {
0470 if (miny[i] < -1E99) {
0471 rh->GetYaxis()->SetRangeUser(miny[i],maxy[i]);
0472 sh->GetYaxis()->SetRangeUser(miny[i],maxy[i]);
0473 }
0474 }
0475
0476
0477
0478
0479 canvas->cd(i+1);
0480
0481 TPad* pad1 = NULL;
0482 TPad* pad2 = NULL;
0483
0484 pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
0485 pad2 = new TPad("pad2", "pad2", 0, 0.0, 1, 0.3);
0486
0487 pad1->SetTopMargin (0.08);
0488 pad1->SetBottomMargin(0.01);
0489 pad1->Draw();
0490
0491 pad2->SetTopMargin (0.05);
0492 pad2->SetBottomMargin(0.45);
0493 pad2->Draw();
0494
0495 pad1->cd();
0496
0497
0498 if (logy) {
0499 if (logy[i])
0500 gPad->SetLogy();
0501
0502 }
0503
0504
0505 #ifdef DEBUG
0506 cout << "DEBUG: Setting statistics..." << endl;
0507 #endif
0508 setStats(sh, rh, -1, 0, false);
0509
0510
0511 #ifdef DEBUG
0512 cout << "DEBUG: Drawing histograms..." << endl;
0513 #endif
0514 if (sh->GetMaximum() > rh->GetMaximum()) {
0515 sh->Draw();
0516 rh->Draw("sames");
0517 }
0518 else {
0519 rh->Draw();
0520 sh->Draw("sames");
0521 }
0522
0523
0524 if (doKolmo) {
0525 #ifdef DEBUG
0526 cout << "DEBUG: Performing Kolmogorov test..." << endl;
0527 #endif
0528 if (doKolmo[i]) {
0529
0530 double kstest = KolmogorovTest(sh,rh);
0531 if(kstest<0.7)
0532 gPad->SetFillColor(kBlue-10);
0533 }
0534 }
0535
0536 pad2->cd();
0537
0538 TH1* ratioplot = PlotRatiosHistograms(rh, sh);
0539 SetHistogramStyle(ratioplot, 21, 4);
0540 ratioplot->Draw("ep");
0541 }
0542
0543
0544
0545 #ifdef DEBUG
0546 cout << "DEBUG: Drawing legend..." << endl;
0547 #endif
0548 canvas->cd();
0549
0550 TLegend* l = 0;
0551 if (nhistos > 4)
0552 l = new TLegend(0.20,0.665,0.80,0.685);
0553 else
0554 l = new TLegend(0.20,0.52,0.80,0.55);
0555
0556 l->SetTextSize(0.011);
0557 l->SetLineColor(1);
0558 l->SetLineWidth(1);
0559 l->SetLineStyle(1);
0560 l->SetFillColor(0);
0561 l->SetBorderSize(2);
0562 l->AddEntry(rh,refLabel,"LPF");
0563 l->AddEntry(sh,newLabel,"LPF");
0564 l->Draw();
0565
0566
0567 canvas->SaveAs(pdfFile+".pdf");
0568 canvas->SaveAs(pdfFile+".png");
0569
0570
0571
0572 delete canvas;
0573 cout << " ... plotted histograms for " << canvasTitle << endl;
0574 }
0575
0576
0577
0578
0579
0580
0581
0582
0583 void Plot4Histograms(const TString& pdfFile,
0584 TDirectory* rdir, TDirectory* sdir,
0585 const TString& rcollname, const TString& scollname,
0586 const char* canvasName, const char* canvasTitle,
0587 const TString& refLabel, const TString& newLabel,
0588 const char** hnames, const char** htitles = 0,
0589 bool* logy = 0, bool* doKolmo = 0, Double_t* norm = 0,bool *resol = 0,
0590 Double_t* minx = 0, Double_t* maxx = 0,
0591 Double_t* miny = 0, Double_t* maxy = 0) {
0592 PlotNHistograms(pdfFile,
0593 rdir, sdir,
0594 rcollname, scollname,
0595 canvasName, canvasTitle,
0596 refLabel, newLabel,
0597 4, hnames, htitles,
0598 logy, doKolmo, norm,resol,minx,maxx,miny,maxy);
0599 }
0600
0601
0602
0603
0604
0605
0606
0607
0608 void Plot6Histograms(const TString& pdfFile,
0609 TDirectory* rdir, TDirectory* sdir,
0610 const TString& rcollname, const TString& scollname,
0611 const char* canvasName, const char* canvasTitle,
0612 const TString& refLabel, const TString& newLabel,
0613 const char** hnames, const char** htitles = 0,
0614 bool* logy = 0, bool* doKolmo = 0, Double_t* norm = 0,bool *resol=0,
0615 Double_t* minx = 0, Double_t* maxx = 0,
0616 Double_t* miny = 0, Double_t* maxy = 0) {
0617
0618
0619 PlotNHistograms(pdfFile,
0620 rdir, sdir,
0621 rcollname, scollname,
0622 canvasName, canvasTitle,
0623 refLabel, newLabel,
0624 6, hnames, htitles,
0625 logy, doKolmo, norm,resol,
0626 minx, maxx,
0627 miny, maxy);
0628 }
0629
0630
0631
0632 void Plot5Histograms(const TString& pdfFile,
0633 TDirectory* rdir, TDirectory* sdir,
0634 const TString& rcollname, const TString& scollname,
0635 const char* canvasName, const char* canvasTitle,
0636 const TString& refLabel, const TString& newLabel,
0637 const char** hnames, const char** htitles = 0,
0638 bool* logy = 0, bool* doKolmo = 0, Double_t* norm = 0,bool *resol=0,
0639 Double_t* minx = 0, Double_t* maxx = 0,
0640 Double_t* miny = 0, Double_t* maxy = 0) {
0641
0642
0643 PlotNHistograms(pdfFile,
0644 rdir, sdir,
0645 rcollname, scollname,
0646 canvasName, canvasTitle,
0647 refLabel, newLabel,
0648 5, hnames, htitles,
0649 logy, doKolmo, norm,resol,
0650 minx, maxx,
0651 miny, maxy);
0652 }
0653
0654
0655
0656
0657
0658
0659
0660
0661 void NormalizeHistogramsToFirst(TH1* h1, TH1* h2) {
0662 if (h1==0 || h2==0) return;
0663
0664 if ( h1->Integral() > 0 && h2->Integral() > 0 ){
0665 Double_t scale2 = h1->Integral()/h2->Integral();
0666 h2->Scale(scale2);
0667 }
0668 }
0669
0670
0671
0672
0673
0674
0675
0676
0677 void NormalizeHistogramsTo1(TH1* h1, TH1* h2) {
0678 if (!h1 || !h2) return;
0679
0680 if ( h1->Integral() != 0 && h2->Integral() != 0 ) {
0681 Double_t scale1 = 1.0/h1->Integral();
0682 Double_t scale2 = 1.0/h2->Integral();
0683 h1->Scale(scale1);
0684 h2->Scale(scale2);
0685 }
0686 }
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
0719
0720
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732
0733
0734
0735
0736
0737
0738
0739
0740
0741
0742
0743
0744
0745
0746
0747
0748
0749
0750
0751
0752
0753
0754
0755
0756
0757
0758
0759
0760
0761
0762
0763
0764
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
0806
0807
0808
0809
0810
0811
0812
0813
0814
0815
0816
0817
0818
0819
0820
0821
0822
0823
0824
0825
0826 TH1* PlotRatiosHistograms(TH1* h1, TH1* h2){
0827
0828 ++ratioCounter;
0829
0830 Int_t nbinsx = h1->GetNbinsX();
0831 Double_t xbins[nbinsx+1];
0832
0833
0834 for (Int_t ibin=1; ibin<=nbinsx+1; ibin++) {
0835 Float_t xmin = h1->GetBinLowEdge(ibin);
0836 xbins[ibin-1] = xmin;
0837 }
0838
0839 TH1F* h_ratio = new TH1F(Form("h_ratio_%d", ratioCounter), "", nbinsx, xbins);
0840
0841 for (Int_t ibin=1; ibin<=nbinsx; ibin++) {
0842
0843 Float_t h1Value = h1->GetBinContent(ibin);
0844 Float_t h1Error = h1->GetBinError(ibin);
0845
0846 Float_t h2Value = h2->GetBinContent(ibin);
0847 Float_t h2Error = h2->GetBinError(ibin);
0848
0849 Float_t ratioVal = 999;
0850 Float_t ratioErr = 999;
0851
0852 if (h2Value > 0 || h2Value < 0) {
0853 ratioVal = h1Value / h2Value;
0854 ratioErr = h1Error / h2Value;
0855 }
0856
0857 h_ratio->SetBinContent(ibin, ratioVal);
0858 h_ratio->SetBinError (ibin, ratioErr);
0859
0860 }
0861
0862 h_ratio->SetTitle("");
0863 h_ratio->GetYaxis()->SetTitle("");
0864 h_ratio->GetYaxis()->SetRangeUser(0.4, 1.6);
0865
0866 return h_ratio;
0867
0868 }