From 7ee672565d4de83c0c8392557767c7e28e55b8d1 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sat, 23 Jul 2022 15:00:24 +0200 Subject: [PATCH 1/2] Update the validation scripts for HFNose in Validation/HGCalValidation and macros in SimG4CMS/Calo --- SimG4CMS/Calo/macros/MakeECPlot.C | 38 ++- SimG4CMS/Calo/macros/MakeGVPlots.C | 44 +--- SimG4CMS/Calo/macros/MakeHFNPlots.C | 22 +- SimG4CMS/Calo/macros/MakeHitStudyPlots.C | 224 +++++++++++++----- .../python/hfnoseDigiStudy_cfi.py | 28 +-- .../python/hfnoseRecHitStudy_cfi.py | 28 +-- .../python/hfnoseSimHitStudy_cfi.py | 27 ++- .../test/python/runHFNoseDigiStudy_cfg.py | 10 +- .../test/python/runHFNoseRecHitStudy_cfg.py | 10 +- .../test/python/runHFNoseSimHitStudy_cfg.py | 8 +- 10 files changed, 267 insertions(+), 172 deletions(-) diff --git a/SimG4CMS/Calo/macros/MakeECPlot.C b/SimG4CMS/Calo/macros/MakeECPlot.C index 282985a3c516b..c1008ed587373 100644 --- a/SimG4CMS/Calo/macros/MakeECPlot.C +++ b/SimG4CMS/Calo/macros/MakeECPlot.C @@ -1,3 +1,33 @@ +////////////////////////////////////////////////////////////////////////////// +// +// Usage: +// .L MakeECPlots.C+g +// +// To make plot from a file created using EcalSimHitStudy +// makeECPlots(fname, text, prefix, save); +// +// where +// fname std::string Name of the ROOT file ("runWithGun_Fix.root") +// text std::string Text written in the Canvas +// ("Fixed Depth Calculation") +// prefix std::string Text added to the name of the Canvas ("Fix") +// save bool Flag to save the canvas as a gif file (false) +// +// To compare plots from 2 different runs with EcalSimHitStudy with +// "old" and "new" settings +// comparePlots(dirname, tex, mom, ratio, fname, save=false) +// +// where +// dirname std::string Name of the directory ("EcalSimHitStudy") +// text std::string Postfix to the histogram name ("All") +// mom int Momentum of the single particle used which +// is also a part of the file name (10) +// ratio bool Shows both distributions or plot the ratio (false) +// fname std:string Prefix of the file name ("elec") +// save bool Flag to save the canvas as a gif file (false) +// +////////////////////////////////////////////////////////////////////////////// + #include #include #include @@ -18,9 +48,9 @@ #include #include -void makePlots(std::string fname="runWithGun_Fix.root", - std::string text="Fixed Depth Calculation", - std::string prefix="Fix", bool save=false) { +void makeECPlots(std::string fname="runWithGun_Fix.root", + std::string text="Fixed Depth Calculation", + std::string prefix="Fix", bool save=false) { std::string name[4] = {"ECLL_EB", "ECLL_EBref", "ECLL_EE", "ECLL_EERef"}; double xrnglo[4] = {1200.0,1200.0,3100.0,3100.0}; double xrnghi[4] = {1600.0,1600.0,3600.0,3600.0}; @@ -49,7 +79,7 @@ void makePlots(std::string fname="runWithGun_Fix.root", } } if (hist != 0) { - sprintf (namep,"%s%s",name[k].c_str(),prefix.c_str()); + sprintf (namep, "%s%s", name[k].c_str(), prefix.c_str()); TCanvas *pad = new TCanvas(namep,namep,500,500); pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); hist->GetYaxis()->SetTitle(ytitl[k].c_str()); diff --git a/SimG4CMS/Calo/macros/MakeGVPlots.C b/SimG4CMS/Calo/macros/MakeGVPlots.C index db79ac614e03c..9d5d61021e037 100644 --- a/SimG4CMS/Calo/macros/MakeGVPlots.C +++ b/SimG4CMS/Calo/macros/MakeGVPlots.C @@ -539,7 +539,6 @@ void setTDRStyle() { // For the Pad: tdrStyle->SetPadBorderMode(0); - // tdrStyle->SetPadBorderSize(Width_t size = 1); tdrStyle->SetPadColor(kWhite); tdrStyle->SetPadGridX(false); tdrStyle->SetPadGridY(false); @@ -557,17 +556,10 @@ void setTDRStyle() { tdrStyle->SetFrameLineWidth(1); // For the histo: - // tdrStyle->SetHistFillColor(1); - // tdrStyle->SetHistFillStyle(0); tdrStyle->SetHistLineColor(1); tdrStyle->SetHistLineStyle(0); tdrStyle->SetHistLineWidth(1); - // tdrStyle->SetLegoInnerR(Float_t rad = 0.5); - // tdrStyle->SetNumberContours(Int_t number = 20); - tdrStyle->SetEndErrorSize(2); - // tdrStyle->SetErrorMarker(20); - //tdrStyle->SetErrorX(0.); tdrStyle->SetMarkerStyle(20); @@ -577,10 +569,9 @@ void setTDRStyle() { tdrStyle->SetFuncColor(2); tdrStyle->SetFuncStyle(1); tdrStyle->SetFuncWidth(1); + //For the date: tdrStyle->SetOptDate(0); - // tdrStyle->SetDateX(Float_t x = 0.01); - // tdrStyle->SetDateY(Float_t y = 0.01); // For the statistics box: tdrStyle->SetOptFile(0); @@ -593,9 +584,6 @@ void setTDRStyle() { tdrStyle->SetStatBorderSize(1); tdrStyle->SetStatH(0.1); tdrStyle->SetStatW(0.15); - // tdrStyle->SetStatStyle(Style_t style = 1001); - // tdrStyle->SetStatX(Float_t x = 0); - // tdrStyle->SetStatY(Float_t y = 0); // Margins: tdrStyle->SetPadTopMargin(0.05); @@ -611,23 +599,14 @@ void setTDRStyle() { tdrStyle->SetTitleTextColor(1); tdrStyle->SetTitleFillColor(10); tdrStyle->SetTitleFontSize(0.05); - // tdrStyle->SetTitleH(0); // Set the height of the title box - // tdrStyle->SetTitleW(0); // Set the width of the title box - // tdrStyle->SetTitleX(0); // Set the position of the title box - // tdrStyle->SetTitleY(0.985); // Set the position of the title box - // tdrStyle->SetTitleStyle(Style_t style = 1001); - // tdrStyle->SetTitleBorderSize(2); // For the axis titles: tdrStyle->SetTitleColor(1, "XYZ"); tdrStyle->SetTitleFont(42, "XYZ"); tdrStyle->SetTitleSize(0.06, "XYZ"); - // tdrStyle->SetTitleXSize(Float_t size = 0.02); // Another way to set the size? - // tdrStyle->SetTitleYSize(Float_t size = 0.02); tdrStyle->SetTitleXOffset(0.9); tdrStyle->SetTitleYOffset(1.25); - // tdrStyle->SetTitleOffset(1.1, "Y"); // Another way to set the Offset // For the axis labels: @@ -652,32 +631,11 @@ void setTDRStyle() { // Postscript options: tdrStyle->SetPaperSize(20.,20.); - // tdrStyle->SetLineScalePS(Float_t scale = 3); - // tdrStyle->SetLineStyleString(Int_t i, const char* text); - // tdrStyle->SetHeaderPS(const char* header); - // tdrStyle->SetTitlePS(const char* pstitle); - - // tdrStyle->SetBarOffset(Float_t baroff = 0.5); - // tdrStyle->SetBarWidth(Float_t barwidth = 0.5); - // tdrStyle->SetPaintTextFormat(const char* format = "g"); - // tdrStyle->SetPalette(Int_t ncolors = 0, Int_t* colors = 0); - // tdrStyle->SetTimeOffset(Double_t toffset); tdrStyle->SetOptLogy(0); tdrStyle->SetOptLogz(0); // Postscript options: tdrStyle->SetPaperSize(20.,20.); - // tdrStyle->SetLineScalePS(Float_t scale = 3); - // tdrStyle->SetLineStyleString(Int_t i, const char* text); - // tdrStyle->SetHeaderPS(const char* header); - // tdrStyle->SetTitlePS(const char* pstitle); - - // tdrStyle->SetBarOffset(Float_t baroff = 0.5); - // tdrStyle->SetBarWidth(Float_t barwidth = 0.5); - // tdrStyle->SetPaintTextFormat(const char* format = "g"); - // tdrStyle->SetPalette(Int_t ncolors = 0, Int_t* colors = 0); - // tdrStyle->SetTimeOffset(Double_t toffset); - // tdrStyle->SetHistMinimumZero(kTRUE); tdrStyle->SetHatchesLineWidth(5); tdrStyle->SetHatchesSpacing(0.05); diff --git a/SimG4CMS/Calo/macros/MakeHFNPlots.C b/SimG4CMS/Calo/macros/MakeHFNPlots.C index e90ddddbabdcd..5f01cec87e989 100644 --- a/SimG4CMS/Calo/macros/MakeHFNPlots.C +++ b/SimG4CMS/Calo/macros/MakeHFNPlots.C @@ -3,19 +3,23 @@ // Usage: // .L MakeHFNPlots.C+g // -// To make plot of layers from a file created using .SimHit/Digi/RecHitStudy +// To make plot of layers from a file created using SimHit/Digi/RecHit +// using the scripts runHFNoseSimHitStudy_cfg.py, runHFNoseDigiStudy_cfg.py, +// runHFNoseRecHitStudy_cfg.py in Validation/HGCalValidation +// // makeLayerPlots(fname, type, todomin, todomax, tag, text, save) // // where -// fname std::string Name of the ROOT file (hfnSimHitD44tt.root) -// type int File type (0:SimHit; 2:Digi; 3:RecHit) -// todomin int Range of plots to make (minimum: -1 -// todomax int maximum: 8) -// tag std::string To be added to the name of the canvas ("") -// text std::string To be added to the title of the histogram ("") -// save bool If the canvas is to be saved as jpg file (false) +// fname std::string Name of the ROOT file [hfnSimHitD94tt.root] +// type int File type (0:SimHit; 2:Digi; 3:RecHit) [0] +// todomin int Range of plots to make [0] +// todomax int (minimum: -1; maximum: 8) [7] +// tag std::string To be added to the name of the canvas [""] +// text std::string To be added to the title of the histogram [""] +// save bool If the canvas is to be saved as jpg file [false] // ////////////////////////////////////////////////////////////////////////////// + #include #include #include @@ -37,7 +41,7 @@ #include #include -void makeLayerPlots(std::string fname="hfnSimHitD44tt.root", int type=0, +void makeLayerPlots(std::string fname="hfnSimHitD94tt.root", int type=0, int todomin=0, int todomax=7, std::string tag="", std::string text="", bool save=false) { diff --git a/SimG4CMS/Calo/macros/MakeHitStudyPlots.C b/SimG4CMS/Calo/macros/MakeHitStudyPlots.C index 7e1dd2546a3cd..1b15785ff5547 100644 --- a/SimG4CMS/Calo/macros/MakeHitStudyPlots.C +++ b/SimG4CMS/Calo/macros/MakeHitStudyPlots.C @@ -1,3 +1,28 @@ +////////////////////////////////////////////////////////////////////////////// +// +// Usage: +// .L MakeHitStudyPlots.C+g +// +// To make plot of various quantities which are created by CaloSimHitStudy +// from one or two different settings +// +// makeHitStudyPlots(file1, tag1, file2, tag2, toDo, ratio, save, dirnm) +// +// where +// file1 std::string Name of the first ROOT file [analRun3Old.root] +// tag1 std::string Tag for the first file [Old] +// file2 std::string Name of the second ROOT file [analRun3New.root] +// tag2 std::string Tag for the second file [New] +// todo int The plot type to be made [0] +// if -1, 6 different types are plotted +// (3, 5, 8, 9, 10, 11) +// ratio bool if the ratio to be plotted [true] +// (works when both files are active) +// save bool If the canvas is to be saved as jpg file [false] +// dirnm std::string Name of the directory [CaloSimHitStudy] +// +////////////////////////////////////////////////////////////////////////////// + #include #include #include @@ -21,9 +46,10 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", std::string tag1 = "Old", std::string file2 = "analRun3New.root", std::string tag2 = "New", - int todo = 0, - std::string dirnm = "CaloSimHitStudy", - bool save = false) { + int toDo = 0, + bool ratio = true, + bool save = false, + std::string dirnm = "CaloSimHitStudy") { std::string names[18] = {"Edep", "EdepEM", "EdepHad", @@ -44,12 +70,17 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", "PtInc"}; int numb[18] = {9, 9, 9, 16, 9, 9, 9, 1, 1, 1, 16, 9, 9, 16, 1, 1, 1, 1}; int rebin[18] = {10, 10, 10, 1, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 1, 1, 1, 1}; + int todos[6] = {3, 5, 8, 9, 10, 11}; + bool debug(false); gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); - gStyle->SetOptStat(1110); + if (ratio) + gStyle->SetOptStat(0); + else + gStyle->SetOptStat(1110); TFile* file[2]; int nfile(0); std::string tag(""), tags[2]; @@ -69,64 +100,141 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", tag += tag2; } } - std::cout << "Use " << nfile << " files from " << file1 << " and " << file2 << std::endl; - for (int i1 = 0; i1 < numb[todo]; ++i1) { - int first(0); - double y1(0.90), dy(0.12); - double y2 = y1 - dy * nfile - 0.01; - TLegend* leg = new TLegend(0.74, y2 - nfile * 0.04, 0.89, y2); - leg->SetBorderSize(1); - leg->SetFillColor(10); - TCanvas* pad; - for (int ifile = 0; ifile < nfile; ++ifile) { - TDirectory* dir = (TDirectory*)file[ifile]->FindObjectAny(dirnm.c_str()); + int todoMin = (toDo >= 0) ? 0 : 0; + int todoMax = (toDo >= 0) ? 0 : 5; + std::cout << "Use " << nfile << " files from " << file1 << " and " << file2 <<" and look for " << todoMin << ":" << todoMax << std::endl; + for (int i0 = todoMin; i0 <= todoMax; ++i0) { + int todo = (todoMax == 0) ? toDo : todos[i0]; + for (int i1 = 0; i1 < numb[todo]; ++i1) { + double y1(0.90), dy(0.12); + double y2 = y1 - dy * nfile - 0.01; + TLegend* leg = (ratio) ? (new TLegend(0.10, 0.86, 0.90, 0.90)) : (new TLegend(0.65, y2 - nfile * 0.04, 0.90, y2)); + leg->SetBorderSize(1); + leg->SetFillColor(10); + TCanvas* pad; + TH1D* hist0[nfile]; char name[100], namec[100]; - if (numb[todo] == 1) { - sprintf(name, "%s", names[todo].c_str()); - sprintf(namec, "%s%s", names[todo].c_str(), tag.c_str()); - } else { - sprintf(name, "%s%d", names[todo].c_str(), i1); - sprintf(namec, "%s%d%s", names[todo].c_str(), i1, tag.c_str()); + for (int ifile = 0; ifile < nfile; ++ifile) { + TDirectory* dir = (TDirectory*)file[ifile]->FindObjectAny(dirnm.c_str()); + if (numb[todo] == 1) { + sprintf(name, "%s", names[todo].c_str()); + sprintf(namec, "%s%s", names[todo].c_str(), tag.c_str()); + } else { + sprintf(name, "%s%d", names[todo].c_str(), i1); + sprintf(namec, "%s%d%s", names[todo].c_str(), i1, tag.c_str()); + } + hist0[ifile] = static_cast(dir->FindObjectAny(name)); + if (debug) + std::cout << name << " read out at " << hist0[ifile] << " for " << tags[ifile] << std::endl; } - TH1D* hist = static_cast(dir->FindObjectAny(name)); - std::cout << name << " read out at " << hist << " for " << tags[ifile] << std::endl; - if (hist != nullptr) { - hist->SetLineColor(first + 1); - hist->SetLineStyle(first + 1); - hist->GetYaxis()->SetTitleOffset(1.4); - if (rebin[todo] > 1) - hist->Rebin(rebin[todo]); - if (first == 0) { - pad = new TCanvas(namec, namec, 500, 500); - pad->SetRightMargin(0.10); - pad->SetTopMargin(0.10); - pad->SetLogy(); - hist->Draw(); - } else { - hist->Draw("sames"); - } - leg->AddEntry(hist, tags[ifile].c_str(), "lp"); - pad->Update(); - ++first; - TPaveStats* st = ((TPaveStats*)hist->GetListOfFunctions()->FindObject("stats")); - if (st != NULL) { - st->SetLineColor(first); - st->SetTextColor(first); - st->SetY1NDC(y1 - dy); - st->SetY2NDC(y1); - st->SetX1NDC(0.65); - st->SetX2NDC(0.90); - y1 -= dy; - } - pad->Modified(); - pad->Update(); - leg->Draw("same"); - pad->Update(); - if (save) { - sprintf(name, "c_%s.pdf", pad->GetName()); - pad->Print(name); + if (!ratio) { + int first(0); + for (int ifile = 0; ifile < nfile; ++ifile) { + TH1D* hist(hist0[ifile]); + if (hist != nullptr) { + hist->SetLineColor(first + 1); + hist->SetLineStyle(first + 1); + hist->GetYaxis()->SetTitleOffset(1.4); + if (rebin[todo] > 1) + hist->Rebin(rebin[todo]); + if (first == 0) { + pad = new TCanvas(namec, namec, 500, 500); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + pad->SetLogy(); + hist->Draw(); + } else { + hist->Draw("sames"); + } + leg->AddEntry(hist, tags[ifile].c_str(), "lp"); + pad->Update(); + ++first; + TPaveStats* st = ((TPaveStats*)hist->GetListOfFunctions()->FindObject("stats")); + if (st != NULL) { + st->SetLineColor(first); + st->SetTextColor(first); + st->SetY1NDC(y1 - dy); + st->SetY2NDC(y1); + st->SetX1NDC(0.65); + st->SetX2NDC(0.90); + y1 -= dy; + } + pad->Modified(); + pad->Update(); + leg->Draw("same"); + pad->Update(); + if (save) { + sprintf(name, "c_%s.pdf", pad->GetName()); + pad->Print(name); + } + } } + } else { + if (nfile == 2) { + int nbin = hist0[0]->GetNbinsX(); + int nbinR = nbin/rebin[todo]; + double xlow = hist0[0]->GetXaxis()->GetBinLowEdge(1); + double xhigh= hist0[0]->GetXaxis()->GetBinLowEdge(nbin) + hist0[0]->GetXaxis()->GetBinWidth(nbin);; + if (numb[todo] == 1) { + sprintf(name, "%sRatio", names[todo].c_str()); + sprintf(namec, "%sRatio%s", names[todo].c_str(), tag.c_str()); + } else { + sprintf(name, "%s%dRatio", names[todo].c_str(), i1); + sprintf(namec, "%s%dRatio%s", names[todo].c_str(), i1, tag.c_str()); + } + pad = new TCanvas(namec, namec, 500, 500); + TH1D* histr = new TH1D(name, hist0[0]->GetTitle(), nbinR, xlow, xhigh); + sprintf (name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str()); + histr->GetXaxis()->SetTitle(hist0[0]->GetXaxis()->GetTitle()); + histr->GetYaxis()->SetTitle(name); + histr->GetXaxis()->SetLabelOffset(0.005); + histr->GetXaxis()->SetTitleOffset(1.00); + histr->GetYaxis()->SetLabelOffset(0.005); + histr->GetYaxis()->SetTitleOffset(1.20); + histr->GetYaxis()->SetRangeUser(0.0,2.0); + double sumNum(0), sumDen(0), maxDev(0); + for (int j=0; jGetBinContent(ib); + tden += hist0[1]->GetBinContent(ib); + rnum += ((hist0[0]->GetBinError(ib))*(hist0[0]->GetBinError(ib))); + rden += ((hist0[1]->GetBinError(ib))*(hist0[1]->GetBinError(ib))); + } + if (tden > 0 && tnum > 0) { + double rat = tnum/tden; + double err = rat*sqrt((rnum/(tnum*tnum))+(rden/(tden*tden))); + histr->SetBinContent(j+1,rat); + histr->SetBinError(j+1,err); + double temp1 = (rat > 1.0) ? 1.0 / rat : rat; + double temp2 = (rat > 1.0) ? err / (rat * rat) : err; + sumNum += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen += (1.0 / (temp2 * temp2)); + if (fabs(1 - temp1) > maxDev) + maxDev = fabs(1 - temp1); + } + } + histr->Draw(); + sprintf (name, "%s vs %s", tag1.c_str(), tag2.c_str()); + leg->AddEntry(histr, name, "lp"); + leg->Draw("same"); + pad->Update(); + TLine* line = new TLine(xlow,1.0,xhigh,1.0); + line->SetLineColor(2); + line->SetLineWidth(2); + line->SetLineStyle(2); + line->Draw("same"); + pad->Modified(); + pad->Update(); + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + if (sumNum == 0) + sumDen = 0; + std::cout << tag1 << " vs " << tag2 << " " << hist0[0]->GetXaxis()->GetTitle() << " Mean deviation " << sumNum << " +- " << sumDen << " maximum " << maxDev << std::endl; + } } } } } + diff --git a/Validation/HGCalValidation/python/hfnoseDigiStudy_cfi.py b/Validation/HGCalValidation/python/hfnoseDigiStudy_cfi.py index bef6bdf327bfc..c005f79ed3aea 100644 --- a/Validation/HGCalValidation/python/hfnoseDigiStudy_cfi.py +++ b/Validation/HGCalValidation/python/hfnoseDigiStudy_cfi.py @@ -3,18 +3,18 @@ from Validation.HGCalValidation.hgcalDigiStudyEE_cfi import * hfnoseDigiStudy = hgcalDigiStudyEE.clone( - detectorName = cms.string("HGCalHFNoseSensitive"), - digiSource = cms.InputTag("hfnoseDigis","HFNose"), - ifNose = cms.untracked.bool(True), - rMin = cms.untracked.double(0), - rMax = cms.untracked.double(150), - zMin = cms.untracked.double(1000), - zMax = cms.untracked.double(1100), - etaMin = cms.untracked.double(2.5), - etaMax = cms.untracked.double(5.5), - nBinR = cms.untracked.int32(150), - nBinZ = cms.untracked.int32(100), - nBinEta = cms.untracked.int32(150), - layers = cms.untracked.int32(8), - ifLayer = cms.untracked.bool(True) + detectorName = "HGCalHFNoseSensitive", + digiSource = "hfnoseDigis : HFNose", + ifNose = True, + rMin = 0, + rMax = 150, + zMin = 1000, + zMax = 1100, + etaMin = 2.5, + etaMax = 5.5, + nBinR = 150, + nBinZ = 100, + nBinEta = 150, + layers = 8, + ifLayer = True ) diff --git a/Validation/HGCalValidation/python/hfnoseRecHitStudy_cfi.py b/Validation/HGCalValidation/python/hfnoseRecHitStudy_cfi.py index c4480570fbba6..ba6dad62394a7 100644 --- a/Validation/HGCalValidation/python/hfnoseRecHitStudy_cfi.py +++ b/Validation/HGCalValidation/python/hfnoseRecHitStudy_cfi.py @@ -3,18 +3,18 @@ from Validation.HGCalValidation.hgcalRecHitStudyEE_cfi import * hfnoseRecHitStudy = hgcalRecHitStudyEE.clone( - detectorName = cms.string("HGCalHFNoseSensitive"), - source = cms.InputTag("HGCalRecHit", "HGCHFNoseRecHits"), - ifNose = cms.untracked.bool(True), - rMin = cms.untracked.double(0), - rMax = cms.untracked.double(150), - zMin = cms.untracked.double(1000), - zMax = cms.untracked.double(1100), - etaMin = cms.untracked.double(2.5), - etaMax = cms.untracked.double(5.5), - nBinR = cms.untracked.int32(150), - nBinZ = cms.untracked.int32(100), - nBinEta = cms.untracked.int32(150), - layers = cms.untracked.int32(8), - ifLayer = cms.untracked.bool(True) + detectorName = "HGCalHFNoseSensitive", + source = "HGCalRecHit : HGCHFNoseRecHits", + ifNose = True, + rMin = 0, + rMax = 150, + zMin = 1000, + zMax = 1100, + etaMin = 2.5, + etaMax = 5.5, + nBinR = 150, + nBinZ = 100, + nBinEta = 150, + layers = 8, + ifLayer = True ) diff --git a/Validation/HGCalValidation/python/hfnoseSimHitStudy_cfi.py b/Validation/HGCalValidation/python/hfnoseSimHitStudy_cfi.py index c912f467ac7e5..d10f128f29089 100644 --- a/Validation/HGCalValidation/python/hfnoseSimHitStudy_cfi.py +++ b/Validation/HGCalValidation/python/hfnoseSimHitStudy_cfi.py @@ -2,16 +2,17 @@ from Validation.HGCalValidation.hgcalSimHitStudy_cfi import * -hgcalSimHitStudy.detectorNames = cms.vstring('HGCalHFNoseSensitive') -hgcalSimHitStudy.caloHitSources = cms.vstring('HFNoseHits') -hgcalSimHitStudy.rMin = cms.untracked.double(0) -hgcalSimHitStudy.rMax = cms.untracked.double(1500) -hgcalSimHitStudy.zMin = cms.untracked.double(10000) -hgcalSimHitStudy.zMax = cms.untracked.double(11000) -hgcalSimHitStudy.etaMin = cms.untracked.double(2.5) -hgcalSimHitStudy.etaMax = cms.untracked.double(5.5) -hgcalSimHitStudy.nBinR = cms.untracked.int32(150) -hgcalSimHitStudy.nBinZ = cms.untracked.int32(100) -hgcalSimHitStudy.nBinEta = cms.untracked.int32(150) -hgcalSimHitStudy.ifNose = cms.untracked.bool(True) -hgcalSimHitStudy.layers = cms.untracked.int32(8) +hgcalSimHitStudy.detectorNames = ['HGCalHFNoseSensitive'] +hgcalSimHitStudy.caloHitSources = ['HFNoseHits'] +hgcalSimHitStudy.rMin = 0 +hgcalSimHitStudy.rMax = 1500 +hgcalSimHitStudy.zMin = 10000 +hgcalSimHitStudy.zMax = 11000 +hgcalSimHitStudy.etaMin = 2.5 +hgcalSimHitStudy.etaMax = 5.5 +hgcalSimHitStudy.nBinR = 50 +hgcalSimHitStudy.nBinZ = 00 +hgcalSimHitStudy.nBinEta = 50 +hgcalSimHitStudy.ifNose = True +hgcalSimHitStudy.layers = 8 +hgcalSimHitStudy.ifLayer = True diff --git a/Validation/HGCalValidation/test/python/runHFNoseDigiStudy_cfg.py b/Validation/HGCalValidation/test/python/runHFNoseDigiStudy_cfg.py index b356caee3dbc9..72579fec557b4 100644 --- a/Validation/HGCalValidation/test/python/runHFNoseDigiStudy_cfg.py +++ b/Validation/HGCalValidation/test/python/runHFNoseDigiStudy_cfg.py @@ -4,8 +4,8 @@ process = cms.Process('PROD',Phase2C11I13M9) process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") -process.load('Configuration.Geometry.GeometryExtended2026D82_cff') -process.load('Configuration.Geometry.GeometryExtended2026DReco_cff') +process.load('Configuration.Geometry.GeometryExtended2026D94_cff') +process.load('Configuration.Geometry.GeometryExtended2026D94Reco_cff') process.load("Configuration.StandardSequences.MagneticField_cff") process.load('FWCore.MessageService.MessageLogger_cfi') process.load('Configuration.StandardSequences.RawToDigi_cff') @@ -19,9 +19,7 @@ process.MessageLogger.HGCalGeom=dict() process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring( - 'file:step2.root', - ) + fileNames = cms.untracked.vstring('file:step2.root') ) process.maxEvents = cms.untracked.PSet( @@ -29,7 +27,7 @@ ) process.TFileService = cms.Service("TFileService", - fileName = cms.string('hfnDigiD82tt.root'), + fileName = cms.string('hfnDigiD94tt.root'), closeFileFast = cms.untracked.bool(True) ) diff --git a/Validation/HGCalValidation/test/python/runHFNoseRecHitStudy_cfg.py b/Validation/HGCalValidation/test/python/runHFNoseRecHitStudy_cfg.py index 497044cd7441b..b6f54a59634e9 100644 --- a/Validation/HGCalValidation/test/python/runHFNoseRecHitStudy_cfg.py +++ b/Validation/HGCalValidation/test/python/runHFNoseRecHitStudy_cfg.py @@ -4,8 +4,8 @@ from Configuration.Eras.Era_Phase2C11I13M9_cff import Phase2C11I13M9 process = cms.Process('PROD',Phase2C11I13M9) -process.load('Configuration.Geometry.GeometryExtended2026D82_cff') -process.load('Configuration.Geometry.GeometryExtended2026D82Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2026D94_cff') +process.load('Configuration.Geometry.GeometryExtended2026D94Reco_cff') process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") process.load('Configuration.StandardSequences.MagneticField_cff') process.load('Configuration.StandardSequences.Services_cff') @@ -20,9 +20,7 @@ process.MessageLogger.HGCalValidation=dict() process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring( - 'file:step3.root', - ) + fileNames = cms.untracked.vstring('file:step3.root') ) process.maxEvents = cms.untracked.PSet( @@ -32,7 +30,7 @@ process.load('Validation.HGCalValidation.hfnoseRecHitStudy_cfi') process.TFileService = cms.Service("TFileService", - fileName = cms.string('hfnRecHitD82tt.root'), + fileName = cms.string('hfnRecHitD94tt.root'), closeFileFast = cms.untracked.bool(True) ) diff --git a/Validation/HGCalValidation/test/python/runHFNoseSimHitStudy_cfg.py b/Validation/HGCalValidation/test/python/runHFNoseSimHitStudy_cfg.py index 99721ffc90cf2..8d7d1541bef1d 100644 --- a/Validation/HGCalValidation/test/python/runHFNoseSimHitStudy_cfg.py +++ b/Validation/HGCalValidation/test/python/runHFNoseSimHitStudy_cfg.py @@ -5,7 +5,7 @@ process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") process.load("IOMC.EventVertexGenerators.VtxSmearedGauss_cfi") -process.load('Configuration.Geometry.GeometryExtended2026D82_cff') +process.load('Configuration.Geometry.GeometryExtended2026D94_cff') process.load("Configuration.StandardSequences.MagneticField_cff") process.load("Configuration.EventContent.EventContent_cff") process.load('FWCore.MessageService.MessageLogger_cfi') @@ -18,9 +18,7 @@ process.MessageLogger.HGCalValidation=dict() process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring( - 'file:step1.root', - ) + fileNames = cms.untracked.vstring('file:step1.root') ) process.maxEvents = cms.untracked.PSet( @@ -28,7 +26,7 @@ ) process.TFileService = cms.Service("TFileService", - fileName = cms.string('hfnSimHitD82tt.root'), + fileName = cms.string('hfnSimHitD94tt.root'), closeFileFast = cms.untracked.bool(True) ) From dec974736caad644ceefdddc8c78e8af395b51c6 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sat, 23 Jul 2022 15:13:51 +0200 Subject: [PATCH 2/2] Code check --- SimG4CMS/Calo/macros/MakeECPlot.C | 375 +++++---- SimG4CMS/Calo/macros/MakeGVPlots.C | 969 +++++++++++++---------- SimG4CMS/Calo/macros/MakeHFNPlots.C | 72 +- SimG4CMS/Calo/macros/MakeHitStudyPlots.C | 238 +++--- 4 files changed, 913 insertions(+), 741 deletions(-) diff --git a/SimG4CMS/Calo/macros/MakeECPlot.C b/SimG4CMS/Calo/macros/MakeECPlot.C index c1008ed587373..145eaf163f4a7 100644 --- a/SimG4CMS/Calo/macros/MakeECPlot.C +++ b/SimG4CMS/Calo/macros/MakeECPlot.C @@ -8,7 +8,7 @@ // // where // fname std::string Name of the ROOT file ("runWithGun_Fix.root") -// text std::string Text written in the Canvas +// text std::string Text written in the Canvas // ("Fixed Depth Calculation") // prefix std::string Text added to the name of the Canvas ("Fix") // save bool Flag to save the canvas as a gif file (false) @@ -48,207 +48,238 @@ #include #include -void makeECPlots(std::string fname="runWithGun_Fix.root", - std::string text="Fixed Depth Calculation", - std::string prefix="Fix", bool save=false) { +void makeECPlots(std::string fname = "runWithGun_Fix.root", + std::string text = "Fixed Depth Calculation", + std::string prefix = "Fix", + bool save = false) { std::string name[4] = {"ECLL_EB", "ECLL_EBref", "ECLL_EE", "ECLL_EERef"}; - double xrnglo[4] = {1200.0,1200.0,3100.0,3100.0}; - double xrnghi[4] = {1600.0,1600.0,3600.0,3600.0}; - std::string xtitl[4] = {"R_{Cyl} (mm)","R_{Cyl} (mm)","|z| (mm)", "|z| (mm)"}; - std::string ytitl[4] = {"# X_{0} (*100)", "# X_{0} (*100)", "# X_{0} (*100)", - "# X_{0} (*100)"}; - std::string title[4] = {"EB (Normal)", "EB (Reflected)", "EE (Normal)", - "EE (Reflected)"}; + double xrnglo[4] = {1200.0, 1200.0, 3100.0, 3100.0}; + double xrnghi[4] = {1600.0, 1600.0, 3600.0, 3600.0}; + std::string xtitl[4] = {"R_{Cyl} (mm)", "R_{Cyl} (mm)", "|z| (mm)", "|z| (mm)"}; + std::string ytitl[4] = {"# X_{0} (*100)", "# X_{0} (*100)", "# X_{0} (*100)", "# X_{0} (*100)"}; + std::string title[4] = {"EB (Normal)", "EB (Reflected)", "EE (Normal)", "EE (Reflected)"}; - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); gStyle->SetOptStat(11110); - TFile *file = new TFile(fname.c_str()); + TFile *file = new TFile(fname.c_str()); if (file) { char namep[100]; - for (int k=0; k<4; ++k) { - TH2D* hist(0); - for (int i=0; i<4; ++i) { - if (i == 0) sprintf (namep, "%s", name[k].c_str()); - else sprintf (namep, "%s;%d", name[k].c_str(),i); - hist = (TH2D*)file->FindObjectAny(name[k].c_str()); - std::cout << namep << " read out at " << hist << std::endl; - if (hist != 0) { - std::cout << "Entries " << hist->GetEntries() << std::endl; - if (hist->GetEntries() > 0) break; - } + for (int k = 0; k < 4; ++k) { + TH2D *hist(0); + for (int i = 0; i < 4; ++i) { + if (i == 0) + sprintf(namep, "%s", name[k].c_str()); + else + sprintf(namep, "%s;%d", name[k].c_str(), i); + hist = (TH2D *)file->FindObjectAny(name[k].c_str()); + std::cout << namep << " read out at " << hist << std::endl; + if (hist != 0) { + std::cout << "Entries " << hist->GetEntries() << std::endl; + if (hist->GetEntries() > 0) + break; + } } if (hist != 0) { - sprintf (namep, "%s%s", name[k].c_str(), prefix.c_str()); - TCanvas *pad = new TCanvas(namep,namep,500,500); - pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); - hist->GetYaxis()->SetTitle(ytitl[k].c_str()); - hist->GetXaxis()->SetTitle(xtitl[k].c_str()); - hist->SetTitle(title[k].c_str()); - hist->GetXaxis()->SetRangeUser(xrnglo[k],xrnghi[k]); - hist->GetYaxis()->SetTitleOffset(1.4); - hist->Draw(); - pad->Update(); - TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats"); - if (st1 != NULL) { - st1->SetY1NDC(0.70); st1->SetY2NDC(0.90); - st1->SetX1NDC(0.65); st1->SetX2NDC(0.90); - } - TPaveText *txt1 = new TPaveText(0.50,0.60,0.90,0.65,"blNDC"); - txt1->SetFillColor(0); - txt1->AddText(text.c_str()); - pad->Update(); - if (save) { - sprintf (namep, "c_%s%s.gif",name[k].c_str(),prefix.c_str()); - pad->Print(namep); - } + sprintf(namep, "%s%s", name[k].c_str(), prefix.c_str()); + TCanvas *pad = new TCanvas(namep, namep, 500, 500); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + hist->GetYaxis()->SetTitle(ytitl[k].c_str()); + hist->GetXaxis()->SetTitle(xtitl[k].c_str()); + hist->SetTitle(title[k].c_str()); + hist->GetXaxis()->SetRangeUser(xrnglo[k], xrnghi[k]); + hist->GetYaxis()->SetTitleOffset(1.4); + hist->Draw(); + pad->Update(); + TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats"); + if (st1 != NULL) { + st1->SetY1NDC(0.70); + st1->SetY2NDC(0.90); + st1->SetX1NDC(0.65); + st1->SetX2NDC(0.90); + } + TPaveText *txt1 = new TPaveText(0.50, 0.60, 0.90, 0.65, "blNDC"); + txt1->SetFillColor(0); + txt1->AddText(text.c_str()); + pad->Update(); + if (save) { + sprintf(namep, "c_%s%s.gif", name[k].c_str(), prefix.c_str()); + pad->Print(namep); + } } } } } -std::vector comparePlots(std::string dirname="EcalSimHitStudy", - std::string text="All", int mom=10, - bool ratio=false, std::string fname="elec", - bool save=false) { +std::vector comparePlots(std::string dirname = "EcalSimHitStudy", + std::string text = "All", + int mom = 10, + bool ratio = false, + std::string fname = "elec", + bool save = false) { + std::vector tcvs; + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); + gStyle->SetOptTitle(0); + gStyle->SetOptFit(0); + if (ratio) + gStyle->SetOptStat(0); + else + gStyle->SetOptStat(1100); - std::vector tcvs; - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); - gStyle->SetOptTitle(0); gStyle->SetOptFit(0); - if (ratio) gStyle->SetOptStat(0); - else gStyle->SetOptStat(1100); - - std::string tags[2] = {"Old", "New"}; - int color[2] = {2,4}; - int marker[2] = {20,21}; - int style[2] = {1,2}; - int rebin[16] = {50,50,50,50, 2, 2, 2, 2, 2, 2,20,20,20,20,20,20}; - int type[16] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; - int edgex[16] = { 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0}; - std::string name1[16] = {"Etot0", "Etot1", "EtotG0", "EtotG1", - "r1by250", "r1by251", "r1by90", "r1by91", - "r9by250", "r9by251", "sEtaEta0", "sEtaEta1", - "sEtaPhi0", "sEtaPhi1", "sPhiPhi0", "sPhiPhi1"}; - char name[100]; - TFile *file[2]; + std::string tags[2] = {"Old", "New"}; + int color[2] = {2, 4}; + int marker[2] = {20, 21}; + int style[2] = {1, 2}; + int rebin[16] = {50, 50, 50, 50, 2, 2, 2, 2, 2, 2, 20, 20, 20, 20, 20, 20}; + int type[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; + int edgex[16] = {0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0}; + std::string name1[16] = {"Etot0", + "Etot1", + "EtotG0", + "EtotG1", + "r1by250", + "r1by251", + "r1by90", + "r1by91", + "r9by250", + "r9by251", + "sEtaEta0", + "sEtaEta1", + "sEtaPhi0", + "sEtaPhi1", + "sPhiPhi0", + "sPhiPhi1"}; + char name[100]; + TFile *file[2]; TDirectory *dir[2]; - for (int i=0; i<2; ++i) { - sprintf (name, "%s%d%s.root", fname.c_str(), mom, tags[i].c_str()); + for (int i = 0; i < 2; ++i) { + sprintf(name, "%s%d%s.root", fname.c_str(), mom, tags[i].c_str()); file[i] = new TFile(name); - dir[i] = (TDirectory*)file[i]->FindObjectAny(dirname.c_str()); + dir[i] = (TDirectory *)file[i]->FindObjectAny(dirname.c_str()); } - for (int k=0; k<16; ++k) { - TH1D* hist[2]; - sprintf (name, "%s", name1[k].c_str()); - for (int i=0; i<2; ++i) { - hist[i] = (TH1D*)dir[i]->FindObjectAny(name); + for (int k = 0; k < 16; ++k) { + TH1D *hist[2]; + sprintf(name, "%s", name1[k].c_str()); + for (int i = 0; i < 2; ++i) { + hist[i] = (TH1D *)dir[i]->FindObjectAny(name); if (hist[i] != 0) { - hist[i]->GetXaxis()->SetLabelOffset(0.005); - hist[i]->GetXaxis()->SetTitleOffset(1.00); - hist[i]->GetYaxis()->SetLabelOffset(0.005); - hist[i]->GetYaxis()->SetTitleOffset(1.20); - hist[i]->SetMarkerStyle(marker[i]); - hist[i]->SetMarkerColor(color[i]); - hist[i]->SetLineColor(color[i]); - hist[i]->SetLineStyle(style[i]); - hist[i]->SetLineWidth(2); + hist[i]->GetXaxis()->SetLabelOffset(0.005); + hist[i]->GetXaxis()->SetTitleOffset(1.00); + hist[i]->GetYaxis()->SetLabelOffset(0.005); + hist[i]->GetYaxis()->SetTitleOffset(1.20); + hist[i]->SetMarkerStyle(marker[i]); + hist[i]->SetMarkerColor(color[i]); + hist[i]->SetLineColor(color[i]); + hist[i]->SetLineStyle(style[i]); + hist[i]->SetLineWidth(2); } } if (hist[0] != 0 && hist[1] != 0) { double ytop(0.90), dy(0.05); double xmin = (edgex[k] == 0) ? 0.65 : 0.11; - double xmin1= (edgex[k] == 0) ? 0.55 : 0.11; - double ymax = ratio ? (ytop - 0.01) : (ytop - 2*dy - 0.01); + double xmin1 = (edgex[k] == 0) ? 0.55 : 0.11; + double ymax = ratio ? (ytop - 0.01) : (ytop - 2 * dy - 0.01); double ymin = ratio ? (ymax - 0.045) : (ymax - 0.09); - TLegend *legend = new TLegend(xmin1, ymin, xmin1+0.35, ymax); + TLegend *legend = new TLegend(xmin1, ymin, xmin1 + 0.35, ymax); legend->SetFillColor(kWhite); if (ratio) { - sprintf (name, "c_R%sE%d%s", name1[k].c_str(), mom, text.c_str()); + sprintf(name, "c_R%sE%d%s", name1[k].c_str(), mom, text.c_str()); } else { - sprintf (name, "c_%sE%d%s", name1[k].c_str(), mom, text.c_str()); + sprintf(name, "c_%sE%d%s", name1[k].c_str(), mom, text.c_str()); } TCanvas *pad = new TCanvas(name, name, 700, 500); pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); - if (type[k] != 0) pad->SetLogy(); + if (type[k] != 0) + pad->SetLogy(); if (ratio) { - int nbin = hist[0]->GetNbinsX(); - int nbinR = nbin/rebin[k]; - double xlow = hist[0]->GetXaxis()->GetBinLowEdge(1); - double xhigh= hist[0]->GetXaxis()->GetBinLowEdge(nbin) + hist[0]->GetXaxis()->GetBinWidth(nbin);; - sprintf (name, "%sRatio", name1[k].c_str()); - TH1D* histr = new TH1D(name, hist[0]->GetTitle(), nbinR, xlow, xhigh); - sprintf (name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str()); - histr->GetXaxis()->SetTitle(hist[0]->GetXaxis()->GetTitle()); - histr->GetYaxis()->SetTitle(name); - histr->GetXaxis()->SetLabelOffset(0.005); - histr->GetXaxis()->SetTitleOffset(1.00); - histr->GetYaxis()->SetLabelOffset(0.005); - histr->GetYaxis()->SetTitleOffset(1.20); - histr->GetYaxis()->SetRangeUser(0.0,2.0); - histr->SetMarkerStyle(marker[0]); - histr->SetMarkerColor(color[0]); - histr->SetLineColor(color[0]); - histr->SetLineStyle(style[0]); - for (int j=0; jGetBinContent(ib); - tden += hist[1]->GetBinContent(ib); - rnum += ((hist[0]->GetBinError(ib))*(hist[0]->GetBinError(ib))); - rden += ((hist[1]->GetBinError(ib))*(hist[1]->GetBinError(ib))); - } - if (tden > 0 && tnum > 0) { - double rat = tnum/tden; - double err = rat*sqrt((rnum/(tnum*tnum))+(rden/(tden*tden))); - histr->SetBinContent(j+1,rat); - histr->SetBinError(j+1,err); - } - } - histr->Draw(); - sprintf (name, "%d GeV Electron (%s)",mom,text.c_str()); - legend->AddEntry(histr,name,"lp"); - pad->Update(); - TLine* line = new TLine(xlow,1.0,xhigh,1.0); - line->SetLineColor(color[1]); - line->SetLineWidth(2); - line->SetLineStyle(2); - line->Draw("same"); - pad->Update(); + int nbin = hist[0]->GetNbinsX(); + int nbinR = nbin / rebin[k]; + double xlow = hist[0]->GetXaxis()->GetBinLowEdge(1); + double xhigh = hist[0]->GetXaxis()->GetBinLowEdge(nbin) + hist[0]->GetXaxis()->GetBinWidth(nbin); + ; + sprintf(name, "%sRatio", name1[k].c_str()); + TH1D *histr = new TH1D(name, hist[0]->GetTitle(), nbinR, xlow, xhigh); + sprintf(name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str()); + histr->GetXaxis()->SetTitle(hist[0]->GetXaxis()->GetTitle()); + histr->GetYaxis()->SetTitle(name); + histr->GetXaxis()->SetLabelOffset(0.005); + histr->GetXaxis()->SetTitleOffset(1.00); + histr->GetYaxis()->SetLabelOffset(0.005); + histr->GetYaxis()->SetTitleOffset(1.20); + histr->GetYaxis()->SetRangeUser(0.0, 2.0); + histr->SetMarkerStyle(marker[0]); + histr->SetMarkerColor(color[0]); + histr->SetLineColor(color[0]); + histr->SetLineStyle(style[0]); + for (int j = 0; j < nbinR; ++j) { + double tnum(0), tden(0), rnum(0), rden(0); + for (int i = 0; i < rebin[k]; ++i) { + int ib = j * rebin[k] + i + 1; + tnum += hist[0]->GetBinContent(ib); + tden += hist[1]->GetBinContent(ib); + rnum += ((hist[0]->GetBinError(ib)) * (hist[0]->GetBinError(ib))); + rden += ((hist[1]->GetBinError(ib)) * (hist[1]->GetBinError(ib))); + } + if (tden > 0 && tnum > 0) { + double rat = tnum / tden; + double err = rat * sqrt((rnum / (tnum * tnum)) + (rden / (tden * tden))); + histr->SetBinContent(j + 1, rat); + histr->SetBinError(j + 1, err); + } + } + histr->Draw(); + sprintf(name, "%d GeV Electron (%s)", mom, text.c_str()); + legend->AddEntry(histr, name, "lp"); + pad->Update(); + TLine *line = new TLine(xlow, 1.0, xhigh, 1.0); + line->SetLineColor(color[1]); + line->SetLineWidth(2); + line->SetLineStyle(2); + line->Draw("same"); + pad->Update(); } else { - double mean[2], error[2]; - for (int i=0; i<2; ++i) { - if (rebin[k] > 1) hist[i]->Rebin(rebin[k]); - if (i == 0) hist[i]->Draw("hist"); - else hist[i]->Draw("sameshist"); - pad->Update(); - sprintf (name, "%d GeV Electron (%s %s)",mom,text.c_str(),tags[i].c_str()); - legend->AddEntry(hist[i],name,"lp"); - TPaveStats* st1 = (TPaveStats*)hist[i]->GetListOfFunctions()->FindObject("stats"); - if (st1 != NULL) { - st1->SetLineColor(color[i]); st1->SetTextColor(color[i]); - st1->SetY1NDC(ytop-dy); st1->SetY2NDC(ytop); - st1->SetX1NDC(xmin); st1->SetX2NDC(xmin+0.25); - ytop -= dy; - } - double entries = hist[i]->GetEntries(); - mean[i] = hist[i]->GetMean(); - error[i] = (hist[i]->GetRMS())/sqrt(entries); - std::cout << text << ":" << hist[i]->GetName() << " V " << tags[i] - << " Mean " << mean[i] << " RMS " << hist[i]->GetRMS() - << " Error " << error[i] << std::endl; - } - double diff = 0.5*(mean[0]-mean[1])/(mean[0]+mean[1]); - double ddiff= (sqrt((mean[0]*error[1])*(mean[0]*error[1])+ - (mean[1]*error[0])*(mean[1]*error[0]))/ - ((mean[0]*mean[0])+(mean[1]*mean[1]))); - double sign = std::abs(diff)/ddiff; - std::cout << "Difference " << diff << " +- " << ddiff - << " Significance " << sign << std::endl; - pad->Modified(); - pad->Update(); + double mean[2], error[2]; + for (int i = 0; i < 2; ++i) { + if (rebin[k] > 1) + hist[i]->Rebin(rebin[k]); + if (i == 0) + hist[i]->Draw("hist"); + else + hist[i]->Draw("sameshist"); + pad->Update(); + sprintf(name, "%d GeV Electron (%s %s)", mom, text.c_str(), tags[i].c_str()); + legend->AddEntry(hist[i], name, "lp"); + TPaveStats *st1 = (TPaveStats *)hist[i]->GetListOfFunctions()->FindObject("stats"); + if (st1 != NULL) { + st1->SetLineColor(color[i]); + st1->SetTextColor(color[i]); + st1->SetY1NDC(ytop - dy); + st1->SetY2NDC(ytop); + st1->SetX1NDC(xmin); + st1->SetX2NDC(xmin + 0.25); + ytop -= dy; + } + double entries = hist[i]->GetEntries(); + mean[i] = hist[i]->GetMean(); + error[i] = (hist[i]->GetRMS()) / sqrt(entries); + std::cout << text << ":" << hist[i]->GetName() << " V " << tags[i] << " Mean " << mean[i] << " RMS " + << hist[i]->GetRMS() << " Error " << error[i] << std::endl; + } + double diff = 0.5 * (mean[0] - mean[1]) / (mean[0] + mean[1]); + double ddiff = + (sqrt((mean[0] * error[1]) * (mean[0] * error[1]) + (mean[1] * error[0]) * (mean[1] * error[0])) / + ((mean[0] * mean[0]) + (mean[1] * mean[1]))); + double sign = std::abs(diff) / ddiff; + std::cout << "Difference " << diff << " +- " << ddiff << " Significance " << sign << std::endl; + pad->Modified(); + pad->Update(); } if (ratio) { } @@ -257,8 +288,8 @@ std::vector comparePlots(std::string dirname="EcalSimHitStudy", pad->Update(); tcvs.push_back(pad); if (save) { - sprintf (name, "%s.pdf", pad->GetName()); - pad->Print(name); + sprintf(name, "%s.pdf", pad->GetName()); + pad->Print(name); } } } diff --git a/SimG4CMS/Calo/macros/MakeGVPlots.C b/SimG4CMS/Calo/macros/MakeGVPlots.C index 9d5d61021e037..7598fb45b1828 100644 --- a/SimG4CMS/Calo/macros/MakeGVPlots.C +++ b/SimG4CMS/Calo/macros/MakeGVPlots.C @@ -8,7 +8,7 @@ // // To plot on the same Canvas plots from Geant4 and GeantV done using // CaloSimHitAnalysis -// makeGV2Plots(fnmG4, fnmGV, todomin, todomax, normalize, tag, text, +// makeGV2Plots(fnmG4, fnmGV, todomin, todomax, normalize, tag, text, // save, modet, mode, dirnm) // // To plot from passive hit collection produced using CaloSimHitAnalysis @@ -31,31 +31,31 @@ // digit indicate text message on the plot // (o for header "CMS Simulation"; t input type; // h for input specification; x for B-field/MT etc) -// mode int Flag if one file is G4 and other GV (0), or +// mode int Flag if one file is G4 and other GV (0), or // both files are GV (1) or G4 (2) (defualt 0) // dirnm std::string Name of the directory ("caloSimHitAnalysis") // // The histogram series have different meanings for each function // // GVPlots (17 in total): -// "Energy deposit per Hit", "Hit time", "Total energy deposit", +// "Energy deposit per Hit", "Hit time", "Total energy deposit", // "Energy deposit after 15 ns", "R vs Z", "R vs Z for hits after 15 ns", // "#phi vs #eta", "Energy deposit", "Time of Hit", -// "Energy deposit per Hit after 15 ns", "Total energy deposit in 100 ns", +// "Energy deposit per Hit after 15 ns", "Total energy deposit in 100 ns", // "Energy deposit for EM particles", "Energy deposit for non-EM particles", // "R", "Z", "#eta", "phi" // // GV2Plots (14 in total): -// "Energy deposit per Hit", "Hit time", "Total energy deposit", +// "Energy deposit per Hit", "Hit time", "Total energy deposit", // "Energy deposit after 15 ns", "Energy deposit", "Time of Hit", -// "Energy deposit per Hit after 15 ns", "Total energy deposit in 100 ns", +// "Energy deposit per Hit after 15 ns", "Total energy deposit in 100 ns", // "Energy deposit for EM particles", "Energy deposit for non-EM particles", // "R", "Z", "#eta", "phi" // // GVSPlots (8 in total): -// "Total Hits", "Tracks", "Step Length', +// "Total Hits", "Tracks", "Step Length', // "Energy Deposit in all components", "Time of all hits", -// "Energy Deposit in tracker subdetectors", +// "Energy Deposit in tracker subdetectors", // "Time of hits in tracker subdetectors" // // All plots in GVPlots or GV2Plots happen for EB, EE, HB and HE @@ -87,454 +87,586 @@ void setTDRStyle(); -void makeGVPlots(std::string fname="analG4.root", bool ifG4=true, - int todomin=0, int todomax=3, std::string tag="", - std::string text="", bool save=false, - std::string dirnm="caloSimHitAnalysis") { - - std::string names[17] = {"EdepT", "TimeT", "Etot", "Edep15", "rz", "rz2", - "etaphi", "Edep", "Time", "EdepT15", "EtotG", - "EdepEM", "EdepHad", "rr", "zz", "eta", "phi"}; - std::string namex[17] = {"Edep", "Time", "Etot", "Edep15", "rz", "rz2", - "etaphi", "EdepX", "TimeX", "EdepT15", "EtotG", - "EdepEM", "EdepHad", "rr", "zz", "eta", "phi"}; - int types[17] = {1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; - int dets[4] = {0, 1, 2, 3}; - std::string detName[4]= {"EB", "EE", "HB", "HE"}; +void makeGVPlots(std::string fname = "analG4.root", + bool ifG4 = true, + int todomin = 0, + int todomax = 3, + std::string tag = "", + std::string text = "", + bool save = false, + std::string dirnm = "caloSimHitAnalysis") { + std::string names[17] = {"EdepT", + "TimeT", + "Etot", + "Edep15", + "rz", + "rz2", + "etaphi", + "Edep", + "Time", + "EdepT15", + "EtotG", + "EdepEM", + "EdepHad", + "rr", + "zz", + "eta", + "phi"}; + std::string namex[17] = {"Edep", + "Time", + "Etot", + "Edep15", + "rz", + "rz2", + "etaphi", + "EdepX", + "TimeX", + "EdepT15", + "EtotG", + "EdepEM", + "EdepHad", + "rr", + "zz", + "eta", + "phi"}; + int types[17] = {1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + int dets[4] = {0, 1, 2, 3}; + std::string detName[4] = {"EB", "EE", "HB", "HE"}; std::string xtitle[17] = {"Energy deposit per Hit (GeV)", - "Hit time (ns)", - "Total energy deposit (GeV)", - "Energy deposit (GeV) after 15 ns", - "z (cm)", "z (cm)", "#eta", - "Energy deposit (GeV)", - "Time of Hit (ns)", - "Energy deposit (GeV) per Hit after 15 ns", - "Total energy deposit (GeV) in 100 ns", - "Energy deposit (GeV) for EM particles", - "Energy deposit (GeV) for non-EM particles", - "R (cm)", "z (cm)", "#eta", "phi"}; - std::string ytitle[17] = {"Hits", "Hits", "Events", "Hits", "R (cm)", - "R (cm)", "#phi", "Hits", "Hits", "Hits", "Events", - "Hits", "Hits", "Hits", "Hits", "Hits", "Hits"}; - - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); + "Hit time (ns)", + "Total energy deposit (GeV)", + "Energy deposit (GeV) after 15 ns", + "z (cm)", + "z (cm)", + "#eta", + "Energy deposit (GeV)", + "Time of Hit (ns)", + "Energy deposit (GeV) per Hit after 15 ns", + "Total energy deposit (GeV) in 100 ns", + "Energy deposit (GeV) for EM particles", + "Energy deposit (GeV) for non-EM particles", + "R (cm)", + "z (cm)", + "#eta", + "phi"}; + std::string ytitle[17] = {"Hits", + "Hits", + "Events", + "Hits", + "R (cm)", + "R (cm)", + "#phi", + "Hits", + "Hits", + "Hits", + "Events", + "Hits", + "Hits", + "Hits", + "Hits", + "Hits", + "Hits"}; + + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); gStyle->SetOptStat(111110); - TFile *file = new TFile(fname.c_str()); + TFile *file = new TFile(fname.c_str()); if (file) { - TDirectory *dir = (TDirectory*)file->FindObjectAny(dirnm.c_str()); + TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnm.c_str()); char cname[100], name[100], name1[100], title[100]; - for (int i1 = 0; i1<4; ++i1) { - for (int i2=todomin; i2<=todomax; ++i2) { - if (types[i2] == 1) { - sprintf (name, "%s%d", names[i2].c_str(), dets[i1]); - sprintf (name1, "%s%d", namex[i2].c_str(), dets[i1]); - if (ifG4) - sprintf (title, "%s %s (Geant4 Simulation)", text.c_str(), detName[i1].c_str()); - else - sprintf (title, "%s %s (GeantV Simulation)", text.c_str(), detName[i1].c_str()); - } else if (i1 == 0) { - sprintf (name, "%s", names[i2].c_str()); - sprintf (name1, "%s", namex[i2].c_str()); - if (ifG4) sprintf (title, "%s (Geant4 Simulation)", text.c_str()); - else sprintf (title, "%s (GeantV Simulation)", text.c_str()); - } else { - continue; - } - TH1D* hist1(nullptr); - TH2D* hist2(nullptr); - if (types[i2] == 1) hist1 = (TH1D*)dir->FindObjectAny(name); - else hist2 = (TH2D*)dir->FindObjectAny(name); -// std::cout << name << " read out at " << hist1 << ":" << hist2 << std::endl; + for (int i1 = 0; i1 < 4; ++i1) { + for (int i2 = todomin; i2 <= todomax; ++i2) { + if (types[i2] == 1) { + sprintf(name, "%s%d", names[i2].c_str(), dets[i1]); + sprintf(name1, "%s%d", namex[i2].c_str(), dets[i1]); + if (ifG4) + sprintf(title, "%s %s (Geant4 Simulation)", text.c_str(), detName[i1].c_str()); + else + sprintf(title, "%s %s (GeantV Simulation)", text.c_str(), detName[i1].c_str()); + } else if (i1 == 0) { + sprintf(name, "%s", names[i2].c_str()); + sprintf(name1, "%s", namex[i2].c_str()); + if (ifG4) + sprintf(title, "%s (Geant4 Simulation)", text.c_str()); + else + sprintf(title, "%s (GeantV Simulation)", text.c_str()); + } else { + continue; + } + TH1D *hist1(nullptr); + TH2D *hist2(nullptr); + if (types[i2] == 1) + hist1 = (TH1D *)dir->FindObjectAny(name); + else + hist2 = (TH2D *)dir->FindObjectAny(name); + // std::cout << name << " read out at " << hist1 << ":" << hist2 << std::endl; if ((hist1 != nullptr) || (hist2 != nullptr)) { - if (ifG4) sprintf (cname, "%sG4%s", name1, tag.c_str()); - else sprintf (cname, "%sGV%s", name1, tag.c_str()); - TCanvas *pad = new TCanvas(cname,cname,500,500); - pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); - if (types[i2] == 1) { - hist1->GetYaxis()->SetTitleOffset(1.2); - hist1->GetYaxis()->SetTitle(ytitle[i2].c_str()); - hist1->GetXaxis()->SetTitle(xtitle[i2].c_str()); - hist1->SetTitle(title); - pad->SetLogy(); - hist1->Draw(); - } else { - hist2->GetYaxis()->SetTitleOffset(1.2); - hist2->GetYaxis()->SetTitle(ytitle[i2].c_str()); - hist2->GetXaxis()->SetTitle(xtitle[i2].c_str()); - hist2->SetTitle(title); - hist2->Draw(); - } - pad->Update(); - TPaveStats* st1 = ((hist1 != nullptr) ? - ((TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats")) : - ((TPaveStats*)hist2->GetListOfFunctions()->FindObject("stats"))); - if (st1 != NULL) { - st1->SetY1NDC(0.70); st1->SetY2NDC(0.90); - st1->SetX1NDC(0.65); st1->SetX2NDC(0.90); - } - pad->Modified(); pad->Update(); - if (save) { - sprintf (name, "c_%s.jpg", pad->GetName()); - pad->Print(name); - } - } + if (ifG4) + sprintf(cname, "%sG4%s", name1, tag.c_str()); + else + sprintf(cname, "%sGV%s", name1, tag.c_str()); + TCanvas *pad = new TCanvas(cname, cname, 500, 500); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + if (types[i2] == 1) { + hist1->GetYaxis()->SetTitleOffset(1.2); + hist1->GetYaxis()->SetTitle(ytitle[i2].c_str()); + hist1->GetXaxis()->SetTitle(xtitle[i2].c_str()); + hist1->SetTitle(title); + pad->SetLogy(); + hist1->Draw(); + } else { + hist2->GetYaxis()->SetTitleOffset(1.2); + hist2->GetYaxis()->SetTitle(ytitle[i2].c_str()); + hist2->GetXaxis()->SetTitle(xtitle[i2].c_str()); + hist2->SetTitle(title); + hist2->Draw(); + } + pad->Update(); + TPaveStats *st1 = ((hist1 != nullptr) ? ((TPaveStats *)hist1->GetListOfFunctions()->FindObject("stats")) + : ((TPaveStats *)hist2->GetListOfFunctions()->FindObject("stats"))); + if (st1 != NULL) { + st1->SetY1NDC(0.70); + st1->SetY2NDC(0.90); + st1->SetX1NDC(0.65); + st1->SetX2NDC(0.90); + } + pad->Modified(); + pad->Update(); + if (save) { + sprintf(name, "c_%s.jpg", pad->GetName()); + pad->Print(name); + } + } } } } } - -void makeGV2Plots(std::string fnmG4="analG4.root", - std::string fnmGV="analGV.root", int todomin=0, - int todomax=2, bool normalize=true, std::string tag="", - std::string text="", int save=0, int modet=0, int mode=0, - std::string dirnm="caloSimHitAnalysis") { - - std::string names[14] = {"EdepT", "TimeT", "Etot", "Edep15", "Edep", "Time", - "EdepT15", "EtotG", "EdepEM", "EdepHad", "rr", "zz", - "eta", "phi"}; - std::string namex[14] = {"Edep", "Time", "Etot", "Edep15", "EdepX", "TimeX", - "EdepT15", "EtotG", "EdepEM", "EdepHad", "rr", "zz", - "eta", "phi"}; - int dets[4] = {0, 1, 2, 3}; +void makeGV2Plots(std::string fnmG4 = "analG4.root", + std::string fnmGV = "analGV.root", + int todomin = 0, + int todomax = 2, + bool normalize = true, + std::string tag = "", + std::string text = "", + int save = 0, + int modet = 0, + int mode = 0, + std::string dirnm = "caloSimHitAnalysis") { + std::string names[14] = {"EdepT", + "TimeT", + "Etot", + "Edep15", + "Edep", + "Time", + "EdepT15", + "EtotG", + "EdepEM", + "EdepHad", + "rr", + "zz", + "eta", + "phi"}; + std::string namex[14] = {"Edep", + "Time", + "Etot", + "Edep15", + "EdepX", + "TimeX", + "EdepT15", + "EtotG", + "EdepEM", + "EdepHad", + "rr", + "zz", + "eta", + "phi"}; + int dets[4] = {0, 1, 2, 3}; std::string detName[4] = {"EB", "EE", "HB", "HE"}; std::string xtitle[14] = {"Energy deposit per Hit (GeV)", - "Hit time (ns)", - "Total energy deposit (GeV)", - "Energy deposit (GeV) after 15 ns", - "Energy deposit (GeV)", - "Time of Hit (ns)", - "Energy deposit (GeV) per Hit after 15 ns", - "Total energy deposit (GeV) in 100 ns", - "Energy deposit (GeV) for EM particles", - "Energy deposit (GeV) for non-EM particles", - "R (cm)", "z (cm)", "#eta", "phi"}; - std::string ytitle[14] = {"Hits", "Hits", "Events", "Hits", "Hits", "Hits", - "Hits", "Events", "Hits", "Hits", "Hits", "Hits", - "Hits", "Hits"}; + "Hit time (ns)", + "Total energy deposit (GeV)", + "Energy deposit (GeV) after 15 ns", + "Energy deposit (GeV)", + "Time of Hit (ns)", + "Energy deposit (GeV) per Hit after 15 ns", + "Total energy deposit (GeV) in 100 ns", + "Energy deposit (GeV) for EM particles", + "Energy deposit (GeV) for non-EM particles", + "R (cm)", + "z (cm)", + "#eta", + "phi"}; + std::string ytitle[14] = {"Hits", + "Hits", + "Events", + "Hits", + "Hits", + "Hits", + "Hits", + "Events", + "Hits", + "Hits", + "Hits", + "Hits", + "Hits", + "Hits"}; std::string title1[4] = {"2 GeV", "10 GeV", "50 GeV", "100 GeV"}; - std::string title2[2] = {"#eta = 1.0; #phi = 1.0", - "|#eta| < 3.0; 0.0 < #phi < 2#pi"}; - std::string title3[4] = {"B = 0.0 Tesla", "B = 3.8 Tesla", - "Multithreaded (B = 0.0 Tesla)", - "Multithreaded (B = 3.8 Tesla)"}; + std::string title2[2] = {"#eta = 1.0; #phi = 1.0", "|#eta| < 3.0; 0.0 < #phi < 2#pi"}; + std::string title3[4] = { + "B = 0.0 Tesla", "B = 3.8 Tesla", "Multithreaded (B = 0.0 Tesla)", "Multithreaded (B = 3.8 Tesla)"}; if (modet >= 10) { setTDRStyle(); } else { - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); } - if (normalize) gStyle->SetOptStat(0); - else gStyle->SetOptStat(111110); - TFile *file1 = new TFile(fnmG4.c_str()); - TFile *file2 = new TFile(fnmGV.c_str()); + if (normalize) + gStyle->SetOptStat(0); + else + gStyle->SetOptStat(111110); + TFile *file1 = new TFile(fnmG4.c_str()); + TFile *file2 = new TFile(fnmGV.c_str()); if (file1 && file2) { - TDirectory *dir1 = (TDirectory*)file1->FindObjectAny(dirnm.c_str()); - TDirectory *dir2 = (TDirectory*)file2->FindObjectAny(dirnm.c_str()); + TDirectory *dir1 = (TDirectory *)file1->FindObjectAny(dirnm.c_str()); + TDirectory *dir2 = (TDirectory *)file2->FindObjectAny(dirnm.c_str()); char name[100], cname[100], title[100]; - for (int i1 = 0; i1<4; ++i1) { - for (int i2=todomin; i2<=todomax; ++i2) { - sprintf (name, "%s%d", names[i2].c_str(), dets[i1]); - sprintf (cname, "%s%s%s", namex[i2].c_str(), detName[i1].c_str(), - tag.c_str()); - if (modet >= 10) - sprintf (title, ""); - else if (mode == 1) - sprintf (title, "%s %s (2 runs of GeantV)", text.c_str(), detName[i1].c_str()); - else if (mode == 2) - sprintf (title, "%s %s (2 runs of Geant4)", text.c_str(), detName[i1].c_str()); - else - sprintf (title, "%s %s (Geant4 vs GeantV)", text.c_str(), detName[i1].c_str()); - TH1D *hist[2]; - hist[0] = (TH1D*)dir1->FindObjectAny(name); - hist[1] = (TH1D*)dir2->FindObjectAny(name); + for (int i1 = 0; i1 < 4; ++i1) { + for (int i2 = todomin; i2 <= todomax; ++i2) { + sprintf(name, "%s%d", names[i2].c_str(), dets[i1]); + sprintf(cname, "%s%s%s", namex[i2].c_str(), detName[i1].c_str(), tag.c_str()); + if (modet >= 10) + sprintf(title, ""); + else if (mode == 1) + sprintf(title, "%s %s (2 runs of GeantV)", text.c_str(), detName[i1].c_str()); + else if (mode == 2) + sprintf(title, "%s %s (2 runs of Geant4)", text.c_str(), detName[i1].c_str()); + else + sprintf(title, "%s %s (Geant4 vs GeantV)", text.c_str(), detName[i1].c_str()); + TH1D *hist[2]; + hist[0] = (TH1D *)dir1->FindObjectAny(name); + hist[1] = (TH1D *)dir2->FindObjectAny(name); if ((hist[0] != nullptr) && (hist[1] != nullptr)) { - // Plot superimposed histograms - double ymx(0.78); - TCanvas *pad = new TCanvas(cname,cname,500,500); - TLegend *legend = new TLegend(0.44, ymx, 0.64, 0.89); - pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); pad->SetLogy(); - pad->SetFillColor(kWhite); legend->SetFillColor(kWhite); - int icol[2] = {1,2}; - int isty[2] = {1,2}; - int imty[2] = {20,24}; - std::string type[2] = {"Geant4", "GeantV"}; - double ymax(0.90); - double total[2] = {0,0}; - for (int i=0; i<2; ++i) { - hist[i]->GetYaxis()->SetTitleOffset(1.2); - hist[i]->GetYaxis()->SetTitle(ytitle[i2].c_str()); - hist[i]->GetXaxis()->SetTitle(xtitle[i2].c_str()); - hist[i]->SetTitle(title); - hist[i]->SetMarkerStyle(imty[i]); - hist[i]->SetMarkerColor(icol[i]); - hist[i]->SetLineColor(icol[i]); - hist[i]->SetLineStyle(isty[i]); - hist[i]->SetNdivisions(505,"X"); - total[i] = hist[i]->GetEntries(); - if (mode == 1) - sprintf (name, "%s (run %d)", type[1].c_str(), i); - else if (mode == 2) - sprintf (name, "%s (run %d)", type[0].c_str(), i); - else - sprintf (name, "%s", type[i].c_str()); - legend->AddEntry(hist[i], name, "lp"); - if (i == 0) { - if (normalize) hist[i]->DrawNormalized("hist"); - else hist[i]->Draw(); - } else { - if (normalize) hist[i]->DrawNormalized("sames hist"); - else hist[i]->Draw("sames"); - } - pad->Update(); - } - legend->Draw("same"); - pad->Modified(); pad->Update(); - for (int i=0; i<2; ++i) { - TPaveStats* st = (TPaveStats*)hist[i]->GetListOfFunctions()->FindObject("stats"); - if (st != NULL) { - st->SetLineColor(icol[i]); st->SetTextColor(icol[i]); - st->SetY1NDC(ymax-0.15); st->SetY2NDC(ymax); - st->SetX1NDC(0.65); st->SetX2NDC(0.90); - ymax -= 0.15; - } - } - pad->Modified(); pad->Update(); - if (ymx > ymax) ymx = ymax; - ymx -= 0.005; - int indx = (modet/10)%10; - double ymi = ymx - 0.05; - if (indx > 0 && indx <= 4) { - sprintf (title, "%s Electrons", title1[indx-1].c_str()); - TPaveText *txt0 = new TPaveText(0.70,ymi,0.895,ymx,"blNDC"); - txt0->SetFillColor(0); - txt0->AddText(title); - txt0->Draw("same"); - ymx -= 0.05; - } - indx = (modet/100)%10; - ymi = ymx - 0.05; - if (indx > 0 && indx <= 2) { - sprintf (title, "%s", title2[indx-1].c_str()); - TPaveText *txt0 = new TPaveText(0.65,ymi,0.895,ymx,"blNDC"); - txt0->SetFillColor(0); - txt0->AddText(title); - txt0->Draw("same"); - ymx -= 0.05; - } - indx = (modet/1000)%10; - ymi = ymx - 0.05; - if (indx > 0 && indx <= 4) { - sprintf (title, "%s", title3[indx-1].c_str()); - TPaveText *txt0 = new TPaveText(0.55,ymi,0.895,ymx,"blNDC"); - txt0->SetFillColor(0); - txt0->AddText(title); - txt0->Draw("same"); - ymx -= 0.05; - } - indx = modet%10; - if (indx > 0) { - sprintf (title, "CMS Simulation"); - TPaveText *txt0 = new TPaveText(0.10,0.91,0.35,0.96,"blNDC"); - txt0->SetFillColor(0); - txt0->AddText(title); - txt0->Draw("same"); - } - if (save != 0) { - if (save > 0) sprintf (name, "c_%s.pdf", pad->GetName()); - else sprintf (name, "c_%s.jpg", pad->GetName()); - pad->Print(name); - } - - // Ratio plot - if (normalize) { - int nbin = hist[0]->GetNbinsX(); - double xmin = hist[0]->GetBinLowEdge(1); - double dx = hist[0]->GetBinWidth(1); - double xmax = xmin + nbin*dx; - double fac = total[1]/total[0]; - int npt = 0; - double sumNum(0), sumDen(0), xpt[200], dxp[200], ypt[200], dyp[200]; - for (int i=0; iGetBinContent(i+1) > 0 && - hist[1]->GetBinContent(i+1) > 0) { - ypt[npt] = (fac*hist[0]->GetBinContent(i+1)/ - hist[1]->GetBinContent(i+1)); - double er1 = hist[0]->GetBinError(i+1)/hist[0]->GetBinContent(i+1); - double er2 = hist[1]->GetBinError(i+1)/hist[1]->GetBinContent(i+1); - dyp[npt] = ypt[npt] * sqrt (er1*er1 + er2*er2); - xpt[npt] = xmin + (i-0.5)*dx; - dxp[npt] = 0; - double temp1 = (ypt[npt]>1.0) ? 1.0/ypt[npt] : ypt[npt]; - double temp2 = (ypt[npt]>1.0) ? dyp[npt]/(ypt[npt]*ypt[npt]) : dyp[npt]; - sumNum += (fabs(1-temp1)/(temp2*temp2)); - sumDen += (1.0/(temp2*temp2)); - ++npt; - } - } - sumNum = (sumDen>0) ? (sumNum/sumDen) : 0; - sumDen = (sumDen>0) ? 1.0/sqrt(sumDen) : 0; - TGraphAsymmErrors *graph = new TGraphAsymmErrors(npt, xpt, ypt, dxp, - dxp, dyp, dyp); - graph->SetMarkerStyle(24); - graph->SetMarkerColor(1); - graph->SetMarkerSize(0.8); - graph->SetLineColor(1); - graph->SetLineWidth(2); - sprintf (name, "%sRatio", pad->GetName()); - TCanvas *canvas = new TCanvas(name, name, 500, 400); - gStyle->SetOptStat(0); gPad->SetTopMargin(0.05); - gPad->SetLeftMargin(0.15); gPad->SetRightMargin(0.025); - gPad->SetBottomMargin(0.20); - TH1F *vFrame = canvas->DrawFrame(xmin, 0.01, xmax, 0.5); - vFrame->GetYaxis()->SetRangeUser(0.4,1.5); - vFrame->GetXaxis()->SetLabelSize(0.035); - vFrame->GetYaxis()->SetLabelSize(0.04); - vFrame->GetXaxis()->SetTitleSize(0.045); - vFrame->GetYaxis()->SetTitleSize(0.045); - vFrame->GetYaxis()->SetTitleOffset(1.2); - vFrame->GetXaxis()->SetRangeUser(xmin, xmax); - vFrame->GetYaxis()->SetTitle("Geant4/GeantV"); - sprintf (name, "%s in %s", xtitle[i2].c_str(), detName[i1].c_str()); - vFrame->GetXaxis()->SetTitle(name); - graph->Draw("P"); - TLine *line = new TLine(xmin, 1.0, xmax, 1.0); - line->SetLineStyle(2); line->SetLineWidth(2); line->SetLineColor(kRed); - line->Draw(); - sprintf (title, "Mean Deviation = %5.3f #pm %5.3f", sumNum, sumDen); - TPaveText* text = new TPaveText(0.16, 0.85, 0.60, 0.90, "brNDC"); - text->AddText(title); text->Draw("same"); - canvas->Modified(); canvas->Update(); - if (save != 0) { - if (save > 0) sprintf (name, "c_%s.pdf", canvas->GetName()); - else sprintf (name, "c_%s.jpg", canvas->GetName()); - canvas->Print(name); - } - } - } + // Plot superimposed histograms + double ymx(0.78); + TCanvas *pad = new TCanvas(cname, cname, 500, 500); + TLegend *legend = new TLegend(0.44, ymx, 0.64, 0.89); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + pad->SetLogy(); + pad->SetFillColor(kWhite); + legend->SetFillColor(kWhite); + int icol[2] = {1, 2}; + int isty[2] = {1, 2}; + int imty[2] = {20, 24}; + std::string type[2] = {"Geant4", "GeantV"}; + double ymax(0.90); + double total[2] = {0, 0}; + for (int i = 0; i < 2; ++i) { + hist[i]->GetYaxis()->SetTitleOffset(1.2); + hist[i]->GetYaxis()->SetTitle(ytitle[i2].c_str()); + hist[i]->GetXaxis()->SetTitle(xtitle[i2].c_str()); + hist[i]->SetTitle(title); + hist[i]->SetMarkerStyle(imty[i]); + hist[i]->SetMarkerColor(icol[i]); + hist[i]->SetLineColor(icol[i]); + hist[i]->SetLineStyle(isty[i]); + hist[i]->SetNdivisions(505, "X"); + total[i] = hist[i]->GetEntries(); + if (mode == 1) + sprintf(name, "%s (run %d)", type[1].c_str(), i); + else if (mode == 2) + sprintf(name, "%s (run %d)", type[0].c_str(), i); + else + sprintf(name, "%s", type[i].c_str()); + legend->AddEntry(hist[i], name, "lp"); + if (i == 0) { + if (normalize) + hist[i]->DrawNormalized("hist"); + else + hist[i]->Draw(); + } else { + if (normalize) + hist[i]->DrawNormalized("sames hist"); + else + hist[i]->Draw("sames"); + } + pad->Update(); + } + legend->Draw("same"); + pad->Modified(); + pad->Update(); + for (int i = 0; i < 2; ++i) { + TPaveStats *st = (TPaveStats *)hist[i]->GetListOfFunctions()->FindObject("stats"); + if (st != NULL) { + st->SetLineColor(icol[i]); + st->SetTextColor(icol[i]); + st->SetY1NDC(ymax - 0.15); + st->SetY2NDC(ymax); + st->SetX1NDC(0.65); + st->SetX2NDC(0.90); + ymax -= 0.15; + } + } + pad->Modified(); + pad->Update(); + if (ymx > ymax) + ymx = ymax; + ymx -= 0.005; + int indx = (modet / 10) % 10; + double ymi = ymx - 0.05; + if (indx > 0 && indx <= 4) { + sprintf(title, "%s Electrons", title1[indx - 1].c_str()); + TPaveText *txt0 = new TPaveText(0.70, ymi, 0.895, ymx, "blNDC"); + txt0->SetFillColor(0); + txt0->AddText(title); + txt0->Draw("same"); + ymx -= 0.05; + } + indx = (modet / 100) % 10; + ymi = ymx - 0.05; + if (indx > 0 && indx <= 2) { + sprintf(title, "%s", title2[indx - 1].c_str()); + TPaveText *txt0 = new TPaveText(0.65, ymi, 0.895, ymx, "blNDC"); + txt0->SetFillColor(0); + txt0->AddText(title); + txt0->Draw("same"); + ymx -= 0.05; + } + indx = (modet / 1000) % 10; + ymi = ymx - 0.05; + if (indx > 0 && indx <= 4) { + sprintf(title, "%s", title3[indx - 1].c_str()); + TPaveText *txt0 = new TPaveText(0.55, ymi, 0.895, ymx, "blNDC"); + txt0->SetFillColor(0); + txt0->AddText(title); + txt0->Draw("same"); + ymx -= 0.05; + } + indx = modet % 10; + if (indx > 0) { + sprintf(title, "CMS Simulation"); + TPaveText *txt0 = new TPaveText(0.10, 0.91, 0.35, 0.96, "blNDC"); + txt0->SetFillColor(0); + txt0->AddText(title); + txt0->Draw("same"); + } + if (save != 0) { + if (save > 0) + sprintf(name, "c_%s.pdf", pad->GetName()); + else + sprintf(name, "c_%s.jpg", pad->GetName()); + pad->Print(name); + } + + // Ratio plot + if (normalize) { + int nbin = hist[0]->GetNbinsX(); + double xmin = hist[0]->GetBinLowEdge(1); + double dx = hist[0]->GetBinWidth(1); + double xmax = xmin + nbin * dx; + double fac = total[1] / total[0]; + int npt = 0; + double sumNum(0), sumDen(0), xpt[200], dxp[200], ypt[200], dyp[200]; + for (int i = 0; i < nbin; ++i) { + if (hist[0]->GetBinContent(i + 1) > 0 && hist[1]->GetBinContent(i + 1) > 0) { + ypt[npt] = (fac * hist[0]->GetBinContent(i + 1) / hist[1]->GetBinContent(i + 1)); + double er1 = hist[0]->GetBinError(i + 1) / hist[0]->GetBinContent(i + 1); + double er2 = hist[1]->GetBinError(i + 1) / hist[1]->GetBinContent(i + 1); + dyp[npt] = ypt[npt] * sqrt(er1 * er1 + er2 * er2); + xpt[npt] = xmin + (i - 0.5) * dx; + dxp[npt] = 0; + double temp1 = (ypt[npt] > 1.0) ? 1.0 / ypt[npt] : ypt[npt]; + double temp2 = (ypt[npt] > 1.0) ? dyp[npt] / (ypt[npt] * ypt[npt]) : dyp[npt]; + sumNum += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen += (1.0 / (temp2 * temp2)); + ++npt; + } + } + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + TGraphAsymmErrors *graph = new TGraphAsymmErrors(npt, xpt, ypt, dxp, dxp, dyp, dyp); + graph->SetMarkerStyle(24); + graph->SetMarkerColor(1); + graph->SetMarkerSize(0.8); + graph->SetLineColor(1); + graph->SetLineWidth(2); + sprintf(name, "%sRatio", pad->GetName()); + TCanvas *canvas = new TCanvas(name, name, 500, 400); + gStyle->SetOptStat(0); + gPad->SetTopMargin(0.05); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.025); + gPad->SetBottomMargin(0.20); + TH1F *vFrame = canvas->DrawFrame(xmin, 0.01, xmax, 0.5); + vFrame->GetYaxis()->SetRangeUser(0.4, 1.5); + vFrame->GetXaxis()->SetLabelSize(0.035); + vFrame->GetYaxis()->SetLabelSize(0.04); + vFrame->GetXaxis()->SetTitleSize(0.045); + vFrame->GetYaxis()->SetTitleSize(0.045); + vFrame->GetYaxis()->SetTitleOffset(1.2); + vFrame->GetXaxis()->SetRangeUser(xmin, xmax); + vFrame->GetYaxis()->SetTitle("Geant4/GeantV"); + sprintf(name, "%s in %s", xtitle[i2].c_str(), detName[i1].c_str()); + vFrame->GetXaxis()->SetTitle(name); + graph->Draw("P"); + TLine *line = new TLine(xmin, 1.0, xmax, 1.0); + line->SetLineStyle(2); + line->SetLineWidth(2); + line->SetLineColor(kRed); + line->Draw(); + sprintf(title, "Mean Deviation = %5.3f #pm %5.3f", sumNum, sumDen); + TPaveText *text = new TPaveText(0.16, 0.85, 0.60, 0.90, "brNDC"); + text->AddText(title); + text->Draw("same"); + canvas->Modified(); + canvas->Update(); + if (save != 0) { + if (save > 0) + sprintf(name, "c_%s.pdf", canvas->GetName()); + else + sprintf(name, "c_%s.jpg", canvas->GetName()); + canvas->Print(name); + } + } + } } } } } - -void makeGVSPlots(std::string fnmG4="analG4.root", - std::string fnmGV="analGV.root", - int todomin=0, int todomax=7, - std::string tag="", std::string text="", bool save=false, - std::string dirnm="caloSimHitAnalysis") { - - std::string names[8] = {"hitp", "trackp", "stepp", "edepp", "timep", "volp", - "edept", "timet"}; - std::string xtitle[8] = {"Hits", "Tracks", "Step Length (cm)", - "Energy Deposit (MeV)", "Time (ns)", - "Volume elements", "Energy Deposit (MeV)", - "Time (ns)"}; - std::string ytitle[8] = {"Events", "Events", "Hits", "Hits", "Hits", - "Events", "Hits", "Hits"}; - std::string detName[6] = {"Barrel Pixel", "Forward Pixel", "TIB", "TID", - "TOB", "TEC"}; +void makeGVSPlots(std::string fnmG4 = "analG4.root", + std::string fnmGV = "analGV.root", + int todomin = 0, + int todomax = 7, + std::string tag = "", + std::string text = "", + bool save = false, + std::string dirnm = "caloSimHitAnalysis") { + std::string names[8] = {"hitp", "trackp", "stepp", "edepp", "timep", "volp", "edept", "timet"}; + std::string xtitle[8] = {"Hits", + "Tracks", + "Step Length (cm)", + "Energy Deposit (MeV)", + "Time (ns)", + "Volume elements", + "Energy Deposit (MeV)", + "Time (ns)"}; + std::string ytitle[8] = {"Events", "Events", "Hits", "Hits", "Hits", "Events", "Hits", "Hits"}; + std::string detName[6] = {"Barrel Pixel", "Forward Pixel", "TIB", "TID", "TOB", "TEC"}; std::string particles[3] = {"Electrons", "Positrons", "Photons"}; int boxp[8] = {0, 1, 0, 0, 0, 0, 0, 0}; int nhis[8] = {-4, -4, -4, 1, 1, 1, 6, 6}; - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); gStyle->SetFillColor(kWhite); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); + gStyle->SetFillColor(kWhite); gStyle->SetOptStat(111110); - TFile *file1 = new TFile(fnmG4.c_str()); - TFile *file2 = new TFile(fnmGV.c_str()); + TFile *file1 = new TFile(fnmG4.c_str()); + TFile *file2 = new TFile(fnmGV.c_str()); if (file1 && file2) { - TDirectory *dir1 = (TDirectory*)file1->FindObjectAny(dirnm.c_str()); - TDirectory *dir2 = (TDirectory*)file2->FindObjectAny(dirnm.c_str()); + TDirectory *dir1 = (TDirectory *)file1->FindObjectAny(dirnm.c_str()); + TDirectory *dir2 = (TDirectory *)file2->FindObjectAny(dirnm.c_str()); char name[100], cname[100], title[100]; - for (int i2=todomin; i2<=todomax; ++i2) { + for (int i2 = todomin; i2 <= todomax; ++i2) { int nhism = (nhis[i2] >= 0) ? nhis[i2] : -nhis[i2]; - for (int i1=0; i1FindObjectAny(name); - hist[1] = (TH1D*)dir2->FindObjectAny(name); - if ((hist[0] != nullptr) && (hist[1] != nullptr)) { - // Plot superimposed histograms - TCanvas *pad = new TCanvas(cname,cname,500,500); - TLegend *legend = new TLegend(0.44, 0.78, 0.64, 0.89); - pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); pad->SetLogy(); - pad->SetFillColor(kWhite); legend->SetFillColor(kWhite); - int icol[2] = {1,2}; - int isty[2] = {1,2}; - int imty[2] = {20,24}; - std::string type[2] = {"Geant4", "GeantV"}; - double ymax(0.90); - double total[2] = {0,0}; - double ymaxv[2] = {0,0}; - for (int i=0; i<2; ++i) { - hist[i]->GetYaxis()->SetTitleOffset(1.2); - hist[i]->GetYaxis()->SetTitle(ytitle[i2].c_str()); - hist[i]->GetXaxis()->SetTitle(xtitle[i2].c_str()); - hist[i]->SetTitle(title); - hist[i]->SetMarkerStyle(imty[i]); - hist[i]->SetMarkerColor(icol[i]); - hist[i]->SetLineColor(icol[i]); - hist[i]->SetLineStyle(isty[i]); - hist[i]->SetNdivisions(505,"X"); - total[i] = hist[i]->GetEntries(); - legend->AddEntry(hist[i],type[i].c_str(),"lp"); - ymaxv[i] = hist[i]->GetMaximum(); - } - int first = (ymaxv[0] > ymaxv[1]) ? 0 : 1; - int next = 1 - first; - hist[first]->Draw(); - hist[next]->Draw("sames"); - pad->Update(); - legend->Draw("same"); - pad->Modified(); pad->Update(); - for (int i=0; i<2; ++i) { - TPaveStats* st = (TPaveStats*)hist[i]->GetListOfFunctions()->FindObject("stats"); - if (st != NULL) { - double xl = (boxp[i2] == 0) ? 0.65 : 0.10; - st->SetLineColor(icol[i]); st->SetTextColor(icol[i]); - st->SetY1NDC(ymax-0.15); st->SetY2NDC(ymax); - st->SetX1NDC(xl); st->SetX2NDC(xl+0.25); - ymax -= 0.15; - } - } - pad->Modified(); pad->Update(); - if (save) { - sprintf (name, "c_%s.jpg", pad->GetName()); - pad->Print(name); - } - } + for (int i1 = 0; i1 < nhism; ++i1) { + if (nhis[i2] <= 1 && i1 == 0) { + sprintf(name, "%s", names[i2].c_str()); + sprintf(cname, "%s%s", names[i2].c_str(), tag.c_str()); + sprintf(title, "%s (Geant4 vs GeantV)", text.c_str()); + } else if (nhis[i2] < 0) { + sprintf(name, "%s%d", names[i2].c_str(), i1 - 1); + sprintf(cname, "%s%d%s", names[i2].c_str(), i1 - 1, tag.c_str()); + sprintf(title, "%s in %s (Geant4 vs GeantV)", particles[i1 - 1].c_str(), text.c_str()); + } else { + sprintf(name, "%s%d", names[i2].c_str(), i1); + sprintf(cname, "%s%d%s", names[i2].c_str(), i1, tag.c_str()); + sprintf(title, "%s in %s (Geant4 vs GeantV)", text.c_str(), detName[i1].c_str()); + } + TH1D *hist[2]; + hist[0] = (TH1D *)dir1->FindObjectAny(name); + hist[1] = (TH1D *)dir2->FindObjectAny(name); + if ((hist[0] != nullptr) && (hist[1] != nullptr)) { + // Plot superimposed histograms + TCanvas *pad = new TCanvas(cname, cname, 500, 500); + TLegend *legend = new TLegend(0.44, 0.78, 0.64, 0.89); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + pad->SetLogy(); + pad->SetFillColor(kWhite); + legend->SetFillColor(kWhite); + int icol[2] = {1, 2}; + int isty[2] = {1, 2}; + int imty[2] = {20, 24}; + std::string type[2] = {"Geant4", "GeantV"}; + double ymax(0.90); + double total[2] = {0, 0}; + double ymaxv[2] = {0, 0}; + for (int i = 0; i < 2; ++i) { + hist[i]->GetYaxis()->SetTitleOffset(1.2); + hist[i]->GetYaxis()->SetTitle(ytitle[i2].c_str()); + hist[i]->GetXaxis()->SetTitle(xtitle[i2].c_str()); + hist[i]->SetTitle(title); + hist[i]->SetMarkerStyle(imty[i]); + hist[i]->SetMarkerColor(icol[i]); + hist[i]->SetLineColor(icol[i]); + hist[i]->SetLineStyle(isty[i]); + hist[i]->SetNdivisions(505, "X"); + total[i] = hist[i]->GetEntries(); + legend->AddEntry(hist[i], type[i].c_str(), "lp"); + ymaxv[i] = hist[i]->GetMaximum(); + } + int first = (ymaxv[0] > ymaxv[1]) ? 0 : 1; + int next = 1 - first; + hist[first]->Draw(); + hist[next]->Draw("sames"); + pad->Update(); + legend->Draw("same"); + pad->Modified(); + pad->Update(); + for (int i = 0; i < 2; ++i) { + TPaveStats *st = (TPaveStats *)hist[i]->GetListOfFunctions()->FindObject("stats"); + if (st != NULL) { + double xl = (boxp[i2] == 0) ? 0.65 : 0.10; + st->SetLineColor(icol[i]); + st->SetTextColor(icol[i]); + st->SetY1NDC(ymax - 0.15); + st->SetY2NDC(ymax); + st->SetX1NDC(xl); + st->SetX2NDC(xl + 0.25); + ymax -= 0.15; + } + } + pad->Modified(); + pad->Update(); + if (save) { + sprintf(name, "c_%s.jpg", pad->GetName()); + pad->Print(name); + } + } } } } } void setTDRStyle() { - TStyle *tdrStyle = new TStyle("tdrStyle","Style for P-TDR"); + TStyle *tdrStyle = new TStyle("tdrStyle", "Style for P-TDR"); // For the canvas: tdrStyle->SetCanvasBorderMode(0); tdrStyle->SetCanvasColor(kWhite); - tdrStyle->SetCanvasDefH(600); //Height of canvas - tdrStyle->SetCanvasDefW(600); //Width of canvas - tdrStyle->SetCanvasDefX(0); //POsition on screen + tdrStyle->SetCanvasDefH(600); //Height of canvas + tdrStyle->SetCanvasDefW(600); //Width of canvas + tdrStyle->SetCanvasDefX(0); //POsition on screen tdrStyle->SetCanvasDefY(0); // For the Pad: @@ -554,15 +686,15 @@ void setTDRStyle() { tdrStyle->SetFrameLineColor(1); tdrStyle->SetFrameLineStyle(1); tdrStyle->SetFrameLineWidth(1); - + // For the histo: tdrStyle->SetHistLineColor(1); tdrStyle->SetHistLineStyle(0); tdrStyle->SetHistLineWidth(1); tdrStyle->SetEndErrorSize(2); - + tdrStyle->SetMarkerStyle(20); - + //For the fit/function: tdrStyle->SetOptFit(1); tdrStyle->SetFitFormat("5.4g"); @@ -575,7 +707,7 @@ void setTDRStyle() { // For the statistics box: tdrStyle->SetOptFile(0); - tdrStyle->SetOptStat(0); // To display the mean and RMS: SetOptStat("mr"); + tdrStyle->SetOptStat(0); // To display the mean and RMS: SetOptStat("mr"); tdrStyle->SetStatColor(kWhite); tdrStyle->SetStatFont(42); tdrStyle->SetStatFontSize(0.025); @@ -630,16 +762,15 @@ void setTDRStyle() { tdrStyle->SetOptLogz(0); // Postscript options: - tdrStyle->SetPaperSize(20.,20.); + tdrStyle->SetPaperSize(20., 20.); tdrStyle->SetOptLogy(0); tdrStyle->SetOptLogz(0); -// Postscript options: - tdrStyle->SetPaperSize(20.,20.); + // Postscript options: + tdrStyle->SetPaperSize(20., 20.); tdrStyle->SetHatchesLineWidth(5); tdrStyle->SetHatchesSpacing(0.05); tdrStyle->cd(); - } diff --git a/SimG4CMS/Calo/macros/MakeHFNPlots.C b/SimG4CMS/Calo/macros/MakeHFNPlots.C index 5f01cec87e989..48ce73d1d4120 100644 --- a/SimG4CMS/Calo/macros/MakeHFNPlots.C +++ b/SimG4CMS/Calo/macros/MakeHFNPlots.C @@ -41,51 +41,59 @@ #include #include -void makeLayerPlots(std::string fname="hfnSimHitD94tt.root", int type=0, - int todomin=0, int todomax=7, std::string tag="", - std::string text="", bool save=false) { - +void makeLayerPlots(std::string fname = "hfnSimHitD94tt.root", + int type = 0, + int todomin = 0, + int todomax = 7, + std::string tag = "", + std::string text = "", + bool save = false) { std::string dirnm[3] = {"hgcalSimHitStudy", "hfnoseDigiStudy", "hfnoseRecHitStudy"}; std::string units[3] = {" (mm)", " (cm)", " (cm)"}; std::string titlx[3] = {"SimHit", "Digi", "RecHit"}; - gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite); - gStyle->SetPadColor(kWhite); + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasColor(kWhite); + gStyle->SetPadColor(kWhite); gStyle->SetOptStat(0); - if (type < 0 || type > 2) type = 0; + if (type < 0 || type > 2) + type = 0; TFile *file = new TFile(fname.c_str()); if (file) { - TDirectory *dir = (TDirectory*)file->FindObjectAny(dirnm[type].c_str()); + TDirectory *dir = (TDirectory *)file->FindObjectAny(dirnm[type].c_str()); char cname[100], name[100], title[100], xtitl[100], ytitl[100]; - for (int i=todomin; i<=todomax; ++i) { + for (int i = todomin; i <= todomax; ++i) { if (i < 0) { - sprintf (name, "RZ_HGCalHFNoseSensitive"); - sprintf (title, "%s (%s)", text.c_str(), titlx[type].c_str()); - sprintf (xtitl, "z %s", units[type].c_str()); - sprintf (ytitl, "R %s", units[type].c_str()); + sprintf(name, "RZ_HGCalHFNoseSensitive"); + sprintf(title, "%s (%s)", text.c_str(), titlx[type].c_str()); + sprintf(xtitl, "z %s", units[type].c_str()); + sprintf(ytitl, "R %s", units[type].c_str()); } else { - sprintf (name, "XY_L%d", i+1); - sprintf (title, "%s (Layer %d %s)", text.c_str(), i+1, titlx[type].c_str()); - sprintf (xtitl, "x %s", units[type].c_str()); - sprintf (ytitl, "y %s", units[type].c_str()); + sprintf(name, "XY_L%d", i + 1); + sprintf(title, "%s (Layer %d %s)", text.c_str(), i + 1, titlx[type].c_str()); + sprintf(xtitl, "x %s", units[type].c_str()); + sprintf(ytitl, "y %s", units[type].c_str()); } - TH2D* hist = (TH2D*)dir->FindObjectAny(name); + TH2D *hist = (TH2D *)dir->FindObjectAny(name); std::cout << name << " read out at " << hist << std::endl; if (hist != nullptr) { - sprintf (cname, "%s%s", name, tag.c_str()); - TCanvas *pad = new TCanvas(cname,cname,500,500); - pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); - hist->GetYaxis()->SetTitleOffset(1.2); - hist->GetYaxis()->SetTitle(ytitl); - hist->GetXaxis()->SetTitle(xtitl); - hist->SetTitle(title); - if (i < 0 && type == 0) hist->GetXaxis()->SetNdivisions(5); - hist->Draw("colz"); - pad->Modified(); pad->Update(); - if (save) { - sprintf (name, "c_%s.jpg", pad->GetName()); - pad->Print(name); - } + sprintf(cname, "%s%s", name, tag.c_str()); + TCanvas *pad = new TCanvas(cname, cname, 500, 500); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + hist->GetYaxis()->SetTitleOffset(1.2); + hist->GetYaxis()->SetTitle(ytitl); + hist->GetXaxis()->SetTitle(xtitl); + hist->SetTitle(title); + if (i < 0 && type == 0) + hist->GetXaxis()->SetNdivisions(5); + hist->Draw("colz"); + pad->Modified(); + pad->Update(); + if (save) { + sprintf(name, "c_%s.jpg", pad->GetName()); + pad->Print(name); + } } } } diff --git a/SimG4CMS/Calo/macros/MakeHitStudyPlots.C b/SimG4CMS/Calo/macros/MakeHitStudyPlots.C index 1b15785ff5547..211d55015e3ca 100644 --- a/SimG4CMS/Calo/macros/MakeHitStudyPlots.C +++ b/SimG4CMS/Calo/macros/MakeHitStudyPlots.C @@ -47,8 +47,8 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", std::string file2 = "analRun3New.root", std::string tag2 = "New", int toDo = 0, - bool ratio = true, - bool save = false, + bool ratio = true, + bool save = false, std::string dirnm = "CaloSimHitStudy") { std::string names[18] = {"Edep", "EdepEM", @@ -102,7 +102,8 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", } int todoMin = (toDo >= 0) ? 0 : 0; int todoMax = (toDo >= 0) ? 0 : 5; - std::cout << "Use " << nfile << " files from " << file1 << " and " << file2 <<" and look for " << todoMin << ":" << todoMax << std::endl; + std::cout << "Use " << nfile << " files from " << file1 << " and " << file2 << " and look for " << todoMin << ":" + << todoMax << std::endl; for (int i0 = todoMin; i0 <= todoMax; ++i0) { int todo = (todoMax == 0) ? toDo : todos[i0]; for (int i1 = 0; i1 < numb[todo]; ++i1) { @@ -115,126 +116,127 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", TH1D* hist0[nfile]; char name[100], namec[100]; for (int ifile = 0; ifile < nfile; ++ifile) { - TDirectory* dir = (TDirectory*)file[ifile]->FindObjectAny(dirnm.c_str()); - if (numb[todo] == 1) { - sprintf(name, "%s", names[todo].c_str()); - sprintf(namec, "%s%s", names[todo].c_str(), tag.c_str()); - } else { - sprintf(name, "%s%d", names[todo].c_str(), i1); - sprintf(namec, "%s%d%s", names[todo].c_str(), i1, tag.c_str()); - } - hist0[ifile] = static_cast(dir->FindObjectAny(name)); - if (debug) - std::cout << name << " read out at " << hist0[ifile] << " for " << tags[ifile] << std::endl; + TDirectory* dir = (TDirectory*)file[ifile]->FindObjectAny(dirnm.c_str()); + if (numb[todo] == 1) { + sprintf(name, "%s", names[todo].c_str()); + sprintf(namec, "%s%s", names[todo].c_str(), tag.c_str()); + } else { + sprintf(name, "%s%d", names[todo].c_str(), i1); + sprintf(namec, "%s%d%s", names[todo].c_str(), i1, tag.c_str()); + } + hist0[ifile] = static_cast(dir->FindObjectAny(name)); + if (debug) + std::cout << name << " read out at " << hist0[ifile] << " for " << tags[ifile] << std::endl; } if (!ratio) { - int first(0); - for (int ifile = 0; ifile < nfile; ++ifile) { - TH1D* hist(hist0[ifile]); - if (hist != nullptr) { - hist->SetLineColor(first + 1); - hist->SetLineStyle(first + 1); - hist->GetYaxis()->SetTitleOffset(1.4); - if (rebin[todo] > 1) - hist->Rebin(rebin[todo]); - if (first == 0) { - pad = new TCanvas(namec, namec, 500, 500); - pad->SetRightMargin(0.10); - pad->SetTopMargin(0.10); - pad->SetLogy(); - hist->Draw(); - } else { - hist->Draw("sames"); - } - leg->AddEntry(hist, tags[ifile].c_str(), "lp"); - pad->Update(); - ++first; - TPaveStats* st = ((TPaveStats*)hist->GetListOfFunctions()->FindObject("stats")); - if (st != NULL) { - st->SetLineColor(first); - st->SetTextColor(first); - st->SetY1NDC(y1 - dy); - st->SetY2NDC(y1); - st->SetX1NDC(0.65); - st->SetX2NDC(0.90); - y1 -= dy; - } - pad->Modified(); - pad->Update(); - leg->Draw("same"); - pad->Update(); - if (save) { - sprintf(name, "c_%s.pdf", pad->GetName()); - pad->Print(name); - } - } + int first(0); + for (int ifile = 0; ifile < nfile; ++ifile) { + TH1D* hist(hist0[ifile]); + if (hist != nullptr) { + hist->SetLineColor(first + 1); + hist->SetLineStyle(first + 1); + hist->GetYaxis()->SetTitleOffset(1.4); + if (rebin[todo] > 1) + hist->Rebin(rebin[todo]); + if (first == 0) { + pad = new TCanvas(namec, namec, 500, 500); + pad->SetRightMargin(0.10); + pad->SetTopMargin(0.10); + pad->SetLogy(); + hist->Draw(); + } else { + hist->Draw("sames"); + } + leg->AddEntry(hist, tags[ifile].c_str(), "lp"); + pad->Update(); + ++first; + TPaveStats* st = ((TPaveStats*)hist->GetListOfFunctions()->FindObject("stats")); + if (st != NULL) { + st->SetLineColor(first); + st->SetTextColor(first); + st->SetY1NDC(y1 - dy); + st->SetY2NDC(y1); + st->SetX1NDC(0.65); + st->SetX2NDC(0.90); + y1 -= dy; + } + pad->Modified(); + pad->Update(); + leg->Draw("same"); + pad->Update(); + if (save) { + sprintf(name, "c_%s.pdf", pad->GetName()); + pad->Print(name); + } + } } } else { - if (nfile == 2) { - int nbin = hist0[0]->GetNbinsX(); - int nbinR = nbin/rebin[todo]; - double xlow = hist0[0]->GetXaxis()->GetBinLowEdge(1); - double xhigh= hist0[0]->GetXaxis()->GetBinLowEdge(nbin) + hist0[0]->GetXaxis()->GetBinWidth(nbin);; - if (numb[todo] == 1) { - sprintf(name, "%sRatio", names[todo].c_str()); - sprintf(namec, "%sRatio%s", names[todo].c_str(), tag.c_str()); - } else { - sprintf(name, "%s%dRatio", names[todo].c_str(), i1); - sprintf(namec, "%s%dRatio%s", names[todo].c_str(), i1, tag.c_str()); - } - pad = new TCanvas(namec, namec, 500, 500); - TH1D* histr = new TH1D(name, hist0[0]->GetTitle(), nbinR, xlow, xhigh); - sprintf (name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str()); - histr->GetXaxis()->SetTitle(hist0[0]->GetXaxis()->GetTitle()); - histr->GetYaxis()->SetTitle(name); - histr->GetXaxis()->SetLabelOffset(0.005); - histr->GetXaxis()->SetTitleOffset(1.00); - histr->GetYaxis()->SetLabelOffset(0.005); - histr->GetYaxis()->SetTitleOffset(1.20); - histr->GetYaxis()->SetRangeUser(0.0,2.0); - double sumNum(0), sumDen(0), maxDev(0); - for (int j=0; jGetBinContent(ib); - tden += hist0[1]->GetBinContent(ib); - rnum += ((hist0[0]->GetBinError(ib))*(hist0[0]->GetBinError(ib))); - rden += ((hist0[1]->GetBinError(ib))*(hist0[1]->GetBinError(ib))); - } - if (tden > 0 && tnum > 0) { - double rat = tnum/tden; - double err = rat*sqrt((rnum/(tnum*tnum))+(rden/(tden*tden))); - histr->SetBinContent(j+1,rat); - histr->SetBinError(j+1,err); - double temp1 = (rat > 1.0) ? 1.0 / rat : rat; - double temp2 = (rat > 1.0) ? err / (rat * rat) : err; - sumNum += (fabs(1 - temp1) / (temp2 * temp2)); - sumDen += (1.0 / (temp2 * temp2)); - if (fabs(1 - temp1) > maxDev) - maxDev = fabs(1 - temp1); - } - } - histr->Draw(); - sprintf (name, "%s vs %s", tag1.c_str(), tag2.c_str()); - leg->AddEntry(histr, name, "lp"); - leg->Draw("same"); - pad->Update(); - TLine* line = new TLine(xlow,1.0,xhigh,1.0); - line->SetLineColor(2); - line->SetLineWidth(2); - line->SetLineStyle(2); - line->Draw("same"); - pad->Modified(); - pad->Update(); - sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; - sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; - if (sumNum == 0) - sumDen = 0; - std::cout << tag1 << " vs " << tag2 << " " << hist0[0]->GetXaxis()->GetTitle() << " Mean deviation " << sumNum << " +- " << sumDen << " maximum " << maxDev << std::endl; - } + if (nfile == 2) { + int nbin = hist0[0]->GetNbinsX(); + int nbinR = nbin / rebin[todo]; + double xlow = hist0[0]->GetXaxis()->GetBinLowEdge(1); + double xhigh = hist0[0]->GetXaxis()->GetBinLowEdge(nbin) + hist0[0]->GetXaxis()->GetBinWidth(nbin); + ; + if (numb[todo] == 1) { + sprintf(name, "%sRatio", names[todo].c_str()); + sprintf(namec, "%sRatio%s", names[todo].c_str(), tag.c_str()); + } else { + sprintf(name, "%s%dRatio", names[todo].c_str(), i1); + sprintf(namec, "%s%dRatio%s", names[todo].c_str(), i1, tag.c_str()); + } + pad = new TCanvas(namec, namec, 500, 500); + TH1D* histr = new TH1D(name, hist0[0]->GetTitle(), nbinR, xlow, xhigh); + sprintf(name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str()); + histr->GetXaxis()->SetTitle(hist0[0]->GetXaxis()->GetTitle()); + histr->GetYaxis()->SetTitle(name); + histr->GetXaxis()->SetLabelOffset(0.005); + histr->GetXaxis()->SetTitleOffset(1.00); + histr->GetYaxis()->SetLabelOffset(0.005); + histr->GetYaxis()->SetTitleOffset(1.20); + histr->GetYaxis()->SetRangeUser(0.0, 2.0); + double sumNum(0), sumDen(0), maxDev(0); + for (int j = 0; j < nbinR; ++j) { + double tnum(0), tden(0), rnum(0), rden(0); + for (int i = 0; i < rebin[todo]; ++i) { + int ib = j * rebin[todo] + i + 1; + tnum += hist0[0]->GetBinContent(ib); + tden += hist0[1]->GetBinContent(ib); + rnum += ((hist0[0]->GetBinError(ib)) * (hist0[0]->GetBinError(ib))); + rden += ((hist0[1]->GetBinError(ib)) * (hist0[1]->GetBinError(ib))); + } + if (tden > 0 && tnum > 0) { + double rat = tnum / tden; + double err = rat * sqrt((rnum / (tnum * tnum)) + (rden / (tden * tden))); + histr->SetBinContent(j + 1, rat); + histr->SetBinError(j + 1, err); + double temp1 = (rat > 1.0) ? 1.0 / rat : rat; + double temp2 = (rat > 1.0) ? err / (rat * rat) : err; + sumNum += (fabs(1 - temp1) / (temp2 * temp2)); + sumDen += (1.0 / (temp2 * temp2)); + if (fabs(1 - temp1) > maxDev) + maxDev = fabs(1 - temp1); + } + } + histr->Draw(); + sprintf(name, "%s vs %s", tag1.c_str(), tag2.c_str()); + leg->AddEntry(histr, name, "lp"); + leg->Draw("same"); + pad->Update(); + TLine* line = new TLine(xlow, 1.0, xhigh, 1.0); + line->SetLineColor(2); + line->SetLineWidth(2); + line->SetLineStyle(2); + line->Draw("same"); + pad->Modified(); + pad->Update(); + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + if (sumNum == 0) + sumDen = 0; + std::cout << tag1 << " vs " << tag2 << " " << hist0[0]->GetXaxis()->GetTitle() << " Mean deviation " << sumNum + << " +- " << sumDen << " maximum " << maxDev << std::endl; + } } } } } -