From f5c009df5e286386c3ab8761638d06940b8d5d83 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Wed, 25 Aug 2021 17:03:45 +0200 Subject: [PATCH 1/2] Update he tools to investigate geometry setups --- .../PrintGeomInfo/test/SimFileCompare.cpp | 29 +++- .../test/python/g4OverlapCheckBigXML_cfg.py | 45 +++++ .../python/g4OverlapCheckBigXMLdd4hep_cfg.py | 47 +++++ Validation/Geometry/macros/MatBudgetVolume.C | 164 +++++++++--------- 4 files changed, 200 insertions(+), 85 deletions(-) create mode 100644 SimG4Core/PrintGeomInfo/test/python/g4OverlapCheckBigXML_cfg.py create mode 100644 SimG4Core/PrintGeomInfo/test/python/g4OverlapCheckBigXMLdd4hep_cfg.py diff --git a/SimG4Core/PrintGeomInfo/test/SimFileCompare.cpp b/SimG4Core/PrintGeomInfo/test/SimFileCompare.cpp index ac8a2af3b0dd8..f0e6932c13936 100644 --- a/SimG4Core/PrintGeomInfo/test/SimFileCompare.cpp +++ b/SimG4Core/PrintGeomInfo/test/SimFileCompare.cpp @@ -275,6 +275,8 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi std::cout << "\nEntries in " << fileFile1 << " and " << fileFile2 << " do not match in the content\n"; const double denmin = 0.0001; int kount1(0), kount2(0); + double difmax1(0), difmax2(0); + std::string nameMax(""); if (type == 0) { const double tol1 = 0.00001; for (auto it1 : matFile1) { @@ -285,6 +287,11 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi 0.5 * (it1.second.radl - it2->second.radl) / std::max(denmin, (it1.second.radl + it2->second.radl)); double idif = 0.5 * (it1.second.intl - it2->second.intl) / std::max(denmin, (it1.second.intl + it2->second.intl)); + if (std::abs(rdif) > difmax1) { + difmax1 = std::abs(rdif); + difmax2 = std::abs(idif); + nameMax = it1.first; + } if ((std::abs(rdif) > tol1) || (std::abs(idif) > tol1)) { ++kount2; std::cout << it1.first << " X0 " << it1.second.radl << ":" << it2->second.radl << ":" << rdif << " #L " @@ -293,7 +300,7 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi } } std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol1 - << " or more\n"; + << " or more; the maximum happens for " << nameMax << " with " << difmax1 << ":" << difmax2 << "\n"; } else if (type == 1) { const double tol2 = 0.0001; for (auto it1 : solidFile1) { @@ -302,6 +309,10 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi ++kount1; double vdif = 0.5 * (it1.second.volume - it2->second.volume) / std::max(denmin, (it1.second.volume + it2->second.volume)); + if (std::abs(vdif) > difmax1) { + difmax1 = std::abs(vdif); + nameMax = it1.first; + } if (std::abs(vdif) > tol2) { ++kount2; std::cout << it1.first << " Volume " << it1.second.volume << ":" << it2->second.volume << ":" << vdif @@ -310,7 +321,7 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi } } std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol2 - << " or more\n"; + << " or more; the maximum happens for " << nameMax << " with " << difmax1 << "\n"; } else if (type == 2) { const double tol3 = 0.0001; for (auto it1 : lvFile1) { @@ -319,6 +330,10 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi ++kount1; double vdif = 0.5 * (it1.second.mass - it2->second.mass) / std::max(denmin, (it1.second.mass + it2->second.mass)); + if (std::abs(vdif) > difmax1) { + difmax1 = std::abs(vdif); + nameMax = it1.first; + } if (std::abs(vdif) > tol3) { ++kount2; std::cout << it1.first << " Mass " << it1.second.mass << ":" << it2->second.mass << ":" << vdif << std::endl; @@ -326,7 +341,7 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi } } std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol3 - << " or more\n"; + << " or more; the maximum happens for " << nameMax << " with " << difmax1 << "\n"; } else { const double tol4 = 0.0001; for (auto it1 : pvFile1) { @@ -336,6 +351,12 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi double xdif = (it1.second.xx - it2->second.xx); double ydif = (it1.second.yy - it2->second.yy); double zdif = (it1.second.zz - it2->second.zz); + double vdif = std::max(std::abs(xdif), std::abs(ydif)); + vdif = std::max(vdif, std::abs(zdif)); + if (vdif > difmax1) { + difmax1 = vdif; + nameMax = it1.first; + } if ((std::abs(xdif) > tol4) || (std::abs(ydif) > tol4) || (std::abs(zdif) > tol4)) { ++kount2; std::cout << it1.first << " x " << it1.second.xx << ":" << it2->second.xx << ":" << xdif << " y " @@ -345,7 +366,7 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi } } std::cout << "\n " << kount2 << " out of " << kount1 << " entries having discrpancies at the level of " << tol4 - << " or more\n"; + << " or more; the maximum happens for " << nameMax << " with " << difmax1 << "\n"; } } diff --git a/SimG4Core/PrintGeomInfo/test/python/g4OverlapCheckBigXML_cfg.py b/SimG4Core/PrintGeomInfo/test/python/g4OverlapCheckBigXML_cfg.py new file mode 100644 index 0000000000000..53233b548fd16 --- /dev/null +++ b/SimG4Core/PrintGeomInfo/test/python/g4OverlapCheckBigXML_cfg.py @@ -0,0 +1,45 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.Eras.Era_Run3_cff import Run3 +process = cms.Process('SIM',Run3) +process.load('SimG4Core.PrintGeomInfo.cmsExtendedGeometry2021_cfi') +process.load("Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cff") +process.load("Geometry.EcalCommonData.ecalSimulationParameters_cff") +process.load("Geometry.HcalCommonData.hcalDDDSimConstants_cff") +process.load("Geometry.HcalCommonData.hcalDDDRecConstants_cfi") +process.load("Geometry.MuonNumbering.muonGeometryConstants_cff") +process.load("Geometry.MuonNumbering.muonOffsetESProducer_cff") + +process.load('FWCore.MessageService.MessageLogger_cfi') + +#if hasattr(process,'MessageLogger'): +# process.MessageLogger.HCalGeom=dict() + +from SimG4Core.PrintGeomInfo.g4TestGeometry_cfi import * +process = checkOverlap(process) + +# enable Geant4 overlap check +process.g4SimHits.CheckGeometry = True + +# Geant4 geometry check +process.g4SimHits.G4CheckOverlap.OutputBaseName = cms.string("cms2021") +process.g4SimHits.G4CheckOverlap.OverlapFlag = cms.bool(True) +process.g4SimHits.G4CheckOverlap.Tolerance = cms.double(0.1) +process.g4SimHits.G4CheckOverlap.Resolution = cms.int32(10000) +process.g4SimHits.G4CheckOverlap.Depth = cms.int32(-1) +# tells if NodeName is G4Region or G4PhysicalVolume +process.g4SimHits.G4CheckOverlap.RegionFlag = cms.bool(False) +# list of names +process.g4SimHits.G4CheckOverlap.NodeNames = cms.vstring('OCMS') +# enable dump gdml file +process.g4SimHits.G4CheckOverlap.gdmlFlag = cms.bool(False) +# if defined a G4PhysicsVolume info is printed +process.g4SimHits.G4CheckOverlap.PVname = '' +# if defined a list of daughter volumes is printed +process.g4SimHits.G4CheckOverlap.LVname = '' + +# extra output files, created if a name is not empty +process.g4SimHits.FileNameField = '' +process.g4SimHits.FileNameGDML = '' +process.g4SimHits.FileNameRegions = '' +# diff --git a/SimG4Core/PrintGeomInfo/test/python/g4OverlapCheckBigXMLdd4hep_cfg.py b/SimG4Core/PrintGeomInfo/test/python/g4OverlapCheckBigXMLdd4hep_cfg.py new file mode 100644 index 0000000000000..2161b487dae3c --- /dev/null +++ b/SimG4Core/PrintGeomInfo/test/python/g4OverlapCheckBigXMLdd4hep_cfg.py @@ -0,0 +1,47 @@ +import FWCore.ParameterSet.Config as cms + +from Configuration.Eras.Era_Run3_dd4hep_cff import Run3_dd4hep +process = cms.Process('G4PrintGeometry',Run3_dd4hep) +process.load("Configuration.Geometry.GeometryDD4hep_cff") +process.load("Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cff") +process.load("Geometry.EcalCommonData.ecalSimulationParameters_cff") +process.load("Geometry.HcalCommonData.hcalDDDSimConstants_cff") +process.load("Geometry.HcalCommonData.hcalDDDRecConstants_cfi") +process.load("Geometry.MuonNumbering.muonGeometryConstants_cff") +process.load("Geometry.MuonNumbering.muonOffsetESProducer_cff") + +process.load('FWCore.MessageService.MessageLogger_cfi') + +#if hasattr(process,'MessageLogger'): +# process.MessageLogger.HCalGeom=dict() + +from SimG4Core.PrintGeomInfo.g4TestGeometry_cfi import * +process = checkOverlap(process) + +# enable Geant4 overlap check +process.g4SimHits.CheckGeometry = True + +process.DDDetectorESProducer.confGeomXMLFiles = cms.FileInPath("SimG4Core/PrintGeomInfo/data/dd4hep/cmsExtendedGeometry2021.xml") + +# Geant4 geometry check +process.g4SimHits.G4CheckOverlap.OutputBaseName = cms.string("cms2021") +process.g4SimHits.G4CheckOverlap.OverlapFlag = cms.bool(True) +process.g4SimHits.G4CheckOverlap.Tolerance = cms.double(0.1) +process.g4SimHits.G4CheckOverlap.Resolution = cms.int32(10000) +process.g4SimHits.G4CheckOverlap.Depth = cms.int32(-1) +# tells if NodeName is G4Region or G4PhysicalVolume +process.g4SimHits.G4CheckOverlap.RegionFlag = cms.bool(False) +# list of names +process.g4SimHits.G4CheckOverlap.NodeNames = cms.vstring('cms:OCMS_1') +# enable dump gdml file +process.g4SimHits.G4CheckOverlap.gdmlFlag = cms.bool(False) +# if defined a G4PhysicsVolume info is printed +process.g4SimHits.G4CheckOverlap.PVname = '' +# if defined a list of daughter volumes is printed +process.g4SimHits.G4CheckOverlap.LVname = '' + +# extra output files, created if a name is not empty +process.g4SimHits.FileNameField = '' +process.g4SimHits.FileNameGDML = '' +process.g4SimHits.FileNameRegions = '' +# diff --git a/Validation/Geometry/macros/MatBudgetVolume.C b/Validation/Geometry/macros/MatBudgetVolume.C index 6e64e0c08c6d5..a793007f20012 100644 --- a/Validation/Geometry/macros/MatBudgetVolume.C +++ b/Validation/Geometry/macros/MatBudgetVolume.C @@ -18,14 +18,14 @@ // Compares material budget plots from 2 different files // // etaPhiPlotComp4(filePreFix, tag, plot, ifEta, debug) -// Compares material budget plots from 4 different files: +// Compares material budget plots from 4 different files: // dddXML, dd4hepXML, dddDB, dd4hepDB // // filePreFix (std::string) Prefix to all 4 file names which will be followed // by one of dddXML/dd4hepXML/dddDB/dd4hepDB strings // and finally with *tag* and ".root" // txt (std::string) Part of the y-title coming after #frac for the plot -// ("{DDD}/{DD4Hep}") +// ("{DDD}{DD4Hep}") // /////////////////////////////////////////////////////////////////////////////// @@ -66,13 +66,13 @@ void etaPhiPlotComp(TString fileName1 = "matbdg_run3.root", std::string plot = "intl", bool ifEta = true, std::string tag = "Run3", - std::string txt = "{DDD}/{DD4Hep}", + std::string txt = "{DDD}/{DD4Hep}", bool debug = false); -void etaPhiPlotComp4(std::string filePreFix = "files/matbdgRun3", - std::string tag = "pre6", - std::string plot = "radl", - bool ifEta = true, - bool debug = false); +void etaPhiPlotComp4(std::string filePreFix = "files/matbdgRun3", + std::string tag = "pre6", + std::string plot = "radl", + bool ifEta = true, + bool debug = false); void setStyle(); const int nlay = 13; @@ -242,8 +242,7 @@ void etaPhi2DPlot(TString fileName, std::string plot, bool drawLeg, double maxEt cc1->Modified(); } -void etaPhiPlotComp( - TString fileName1, TString fileName2, std::string plot, bool ifEta, std::string tag, std::string txt, bool debug) { +void etaPhiPlotComp(TString fileName1, TString fileName2, std::string plot, bool ifEta, std::string tag, std::string txt, bool debug) { setStyle(); gStyle->SetOptTitle(0); TFile *file1 = new TFile(fileName1); @@ -261,12 +260,12 @@ void etaPhiPlotComp( std::string ztit = "Eta"; char ytit[40]; if (plot == "radl") { - sprintf(ytit, "#frac%s for MB (X_{0})", txt.c_str()); + sprintf (ytit, "#frac%s for MB (X_{0})", txt.c_str()); } else if (plot == "step") { - sprintf(ytit, "#frac%s for MB (Step Length)", txt.c_str()); + sprintf (ytit, "#frac%s for MB (Step Length)", txt.c_str()); } else { plot = "intl"; - sprintf(ytit, "#frac%s for MB (#lambda)", txt.c_str()); + sprintf (ytit, "#frac%s for MB (#lambda)", txt.c_str()); } if (!ifEta) { xtit = "#phi"; @@ -299,7 +298,7 @@ void etaPhiPlotComp( } std::vector xx, yy, dx, dy; int ii = nflayer[i]; - double sumNum(0), sumDen(0); + double sumNum(0), sumDen(0), maxtmp(0), maxDev(0), dmaxDev(0); for (unsigned int k = 0; k < xx0.size(); ++k) { if ((yy1[k] > 0) && (yy2[k] > 0)) { double rat = yy1[k] / yy2[k]; @@ -314,13 +313,19 @@ void etaPhiPlotComp( } double temp1 = (rat > 1.0) ? 1.0 / rat : rat; double temp2 = (rat > 1.0) ? drt / (rat * rat) : drt; - sumNum += (fabs(1.0 - temp1) / (temp2 * temp2)); + double temp0 = (fabs(1.0 - temp1) / (temp2 * temp2)); + sumNum += temp0; sumDen += (1.0 / (temp2 * temp2)); + if (temp0 >= maxtmp) { + maxtmp = temp0; + maxDev = fabs(1.0 - temp1); + dmaxDev = temp2; + } } } sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; - std::cout << "Mean deviation for " << title[ii] << " " << sumNum << " +- " << sumDen << std::endl; + std::cout << "Mean deviation for " << title[ii] << " " << sumNum << " +- " << sumDen << " Max " << maxDev << " +- " << dmaxDev << std::endl; if (xx.size() > 0) { TGraphErrors *graph = new TGraphErrors(xx.size(), &xx[0], &yy[0], &dx[0], &dy[0]); graph->SetLineColor(colorLay[ii]); @@ -416,74 +421,71 @@ void etaPhiPlotComp4(std::string filePreFix, std::string tag, std::string plot, int ii = nflayer[i] + j; sprintf(hname, "%s%s%s", plot.c_str(), ztit.c_str(), names[ii].c_str()); TProfile *prof[files]; - bool okf(true); - for (int k1 = 0; k1 < files; ++k1) { - dir[k1]->GetObject(hname, prof[k1]); - if (dir[k1] == nullptr) - okf = false; - } + bool okf(true); + for (int k1 = 0; k1 < files; ++k1) { + dir[k1]->GetObject(hname, prof[k1]); + if (dir[k1] == nullptr) + okf = false; + } if (okf) { int nb = prof[0]->GetNbinsX(); for (int k = 1; k <= nb; ++k) { xx0.push_back(prof[0]->GetBinLowEdge(k) + prof[0]->GetBinWidth(k)); - for (int k1 = 0; k1 < files; ++k1) { - yy0[k1].push_back(prof[k1]->GetBinContent(k)); - dy0[k1].push_back(prof[k1]->GetBinError(k)); - } + for (int k1 = 0; k1 < files; ++k1) { + yy0[k1].push_back(prof[k1]->GetBinContent(k)); + dy0[k1].push_back(prof[k1]->GetBinError(k)); + } } } } int ii = nflayer[i]; for (int k1 = 1; k1 < files; ++k1) { - std::vector xx, yy, dx, dy; - double sumNum(0), sumDen(0), maxtmp(0), maxDev(0), dmaxDev(0); - for (unsigned int k = 0; k < xx0.size(); ++k) { - if ((yy0[0][k] > 0) && (yy0[k1][k] > 0)) { - double rat = yy0[k1][k] / yy0[0][k]; - double drt = rat * sqrt((dy0[k1][k] / yy0[k1][k]) * (dy0[k1][k] / yy0[k1][k]) + - (dy0[0][k] / yy0[0][k]) * (dy0[0][k] / yy0[0][k])); - xx.push_back(xx0[k]); - dx.push_back(0); - yy.push_back(rat); - dy.push_back(drt); - if (debug) { - std::cout << nametype[k1] << ":" << title[ii] << " [" << (xx.size() - 1) << "] " << xx0[k] << " Ratio " - << rat << " +- " << drt << std::endl; - } - double temp1 = (rat > 1.0) ? 1.0 / rat : rat; - double temp2 = (rat > 1.0) ? drt / (rat * rat) : drt; - double temp0 = (fabs(1.0 - temp1) / (temp2 * temp2)); - sumNum += temp0; - sumDen += (1.0 / (temp2 * temp2)); - if (temp0 >= maxtmp) { - maxtmp = temp0; - maxDev = fabs(1.0 - temp1); - dmaxDev = temp2; - } - } - } - sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; - sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; - std::cout << title[ii] << " in " << nametype[k1] << " Mean " << sumNum << " +- " << sumDen << " Max " << maxDev - << " +- " << dmaxDev << std::endl; - if (xx.size() > 0) { - TGraphErrors *graph = new TGraphErrors(xx.size(), &xx[0], &yy[0], &dx[0], &dy[0]); - graph->SetLineColor(colortype[k1]); - graph->SetFillColor(colorLay[ii]); - graph->SetMarkerStyle(styleLay[ii]); - if (k1 == 1) { - sprintf(titlex, "%s", title[ii].c_str()); - leg->AddEntry(graph, titlex, "lep"); - } - graphs.push_back(graph); - if (nb == 0) { - sprintf(hname, "%s%s%s", plot.c_str(), ztit.c_str(), names[0].c_str()); - TProfile *prof; - dir[0]->GetObject(hname, prof); - nb = prof->GetNbinsX(); - xlow = prof->GetBinLowEdge(1); - xhigh = prof->GetBinLowEdge(nb) + prof->GetBinWidth(nb); + std::vector xx, yy, dx, dy; + double sumNum(0), sumDen(0), maxtmp(0), maxDev(0), dmaxDev(0); + for (unsigned int k = 0; k < xx0.size(); ++k) { + if ((yy0[0][k] > 0) && (yy0[k1][k] > 0)) { + double rat = yy0[k1][k] / yy0[0][k]; + double drt = rat * sqrt((dy0[k1][k] / yy0[k1][k]) * (dy0[k1][k] / yy0[k1][k]) + (dy0[0][k] / yy0[0][k]) * (dy0[0][k] / yy0[0][k])); + xx.push_back(xx0[k]); + dx.push_back(0); + yy.push_back(rat); + dy.push_back(drt); + if (debug) { + std::cout << nametype[k1] << ":" << title[ii] << " [" << (xx.size() - 1) << "] " << xx0[k] << " Ratio " << rat << " +- " << drt << std::endl; } + double temp1 = (rat > 1.0) ? 1.0 / rat : rat; + double temp2 = (rat > 1.0) ? drt / (rat * rat) : drt; + double temp0 = (fabs(1.0 - temp1) / (temp2 * temp2)); + sumNum += temp0; + sumDen += (1.0 / (temp2 * temp2)); + if (temp0 >= maxtmp) { + maxtmp = temp0; + maxDev = fabs(1.0 - temp1); + dmaxDev = temp2; + } + } + } + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + std::cout << title[ii] << " in " << nametype[k1] << " Mean " << sumNum << " +- " << sumDen << " Max " << maxDev << " +- " << dmaxDev << std::endl; + if (xx.size() > 0) { + TGraphErrors *graph = new TGraphErrors(xx.size(), &xx[0], &yy[0], &dx[0], &dy[0]); + graph->SetLineColor(colortype[k1]); + graph->SetFillColor(colorLay[ii]); + graph->SetMarkerStyle(styleLay[ii]); + if (k1 == 1) { + sprintf(titlex, "%s", title[ii].c_str()); + leg->AddEntry(graph, titlex, "lep"); + } + graphs.push_back(graph); + if (nb == 0) { + sprintf(hname, "%s%s%s", plot.c_str(), ztit.c_str(), names[0].c_str()); + TProfile *prof; + dir[0]->GetObject(hname, prof); + nb = prof->GetNbinsX(); + xlow = prof->GetBinLowEdge(1); + xhigh = prof->GetBinLowEdge(nb) + prof->GetBinWidth(nb); + } } } } @@ -509,13 +511,13 @@ void etaPhiPlotComp4(std::string filePreFix, std::string tag, std::string plot, cc1->Modified(); double ymx = 0.68; for (int k1 = 1; k1 < files; ++k1) { - TPaveText *txt1 = new TPaveText(0.84, ymx - 0.03, 0.99, ymx, "blNDC"); - txt1->SetFillColor(0); - sprintf(fname, "%s", nametype[k1].c_str()); - txt1->AddText(fname); - ((TText *)txt1->GetListOfLines()->Last())->SetTextColor(colortype[k1]); - txt1->Draw(); - ymx -= 0.03; + TPaveText* txt1 = new TPaveText(0.84, ymx - 0.03, 0.99, ymx, "blNDC"); + txt1->SetFillColor(0); + sprintf(fname, "%s", nametype[k1].c_str()); + txt1->AddText(fname); + ((TText*)txt1->GetListOfLines()->Last())->SetTextColor(colortype[k1]); + txt1->Draw(); + ymx -= 0.03; } } } From 8f01892b904ee56f45fbcffec46bf09f3a535f1c Mon Sep 17 00:00:00 2001 From: Sunanda Date: Wed, 25 Aug 2021 17:18:46 +0200 Subject: [PATCH 2/2] Code check --- .../PrintGeomInfo/test/SimFileCompare.cpp | 38 ++-- Validation/Geometry/macros/MatBudgetVolume.C | 165 +++++++++--------- 2 files changed, 104 insertions(+), 99 deletions(-) diff --git a/SimG4Core/PrintGeomInfo/test/SimFileCompare.cpp b/SimG4Core/PrintGeomInfo/test/SimFileCompare.cpp index f0e6932c13936..bea0d0746eec2 100644 --- a/SimG4Core/PrintGeomInfo/test/SimFileCompare.cpp +++ b/SimG4Core/PrintGeomInfo/test/SimFileCompare.cpp @@ -287,11 +287,11 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi 0.5 * (it1.second.radl - it2->second.radl) / std::max(denmin, (it1.second.radl + it2->second.radl)); double idif = 0.5 * (it1.second.intl - it2->second.intl) / std::max(denmin, (it1.second.intl + it2->second.intl)); - if (std::abs(rdif) > difmax1) { - difmax1 = std::abs(rdif); - difmax2 = std::abs(idif); - nameMax = it1.first; - } + if (std::abs(rdif) > difmax1) { + difmax1 = std::abs(rdif); + difmax2 = std::abs(idif); + nameMax = it1.first; + } if ((std::abs(rdif) > tol1) || (std::abs(idif) > tol1)) { ++kount2; std::cout << it1.first << " X0 " << it1.second.radl << ":" << it2->second.radl << ":" << rdif << " #L " @@ -309,10 +309,10 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi ++kount1; double vdif = 0.5 * (it1.second.volume - it2->second.volume) / std::max(denmin, (it1.second.volume + it2->second.volume)); - if (std::abs(vdif) > difmax1) { - difmax1 = std::abs(vdif); - nameMax = it1.first; - } + if (std::abs(vdif) > difmax1) { + difmax1 = std::abs(vdif); + nameMax = it1.first; + } if (std::abs(vdif) > tol2) { ++kount2; std::cout << it1.first << " Volume " << it1.second.volume << ":" << it2->second.volume << ":" << vdif @@ -330,10 +330,10 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi ++kount1; double vdif = 0.5 * (it1.second.mass - it2->second.mass) / std::max(denmin, (it1.second.mass + it2->second.mass)); - if (std::abs(vdif) > difmax1) { - difmax1 = std::abs(vdif); - nameMax = it1.first; - } + if (std::abs(vdif) > difmax1) { + difmax1 = std::abs(vdif); + nameMax = it1.first; + } if (std::abs(vdif) > tol3) { ++kount2; std::cout << it1.first << " Mass " << it1.second.mass << ":" << it2->second.mass << ":" << vdif << std::endl; @@ -351,12 +351,12 @@ void CompareFiles(const char* fileFile1, const char* fileFile2, int type, int fi double xdif = (it1.second.xx - it2->second.xx); double ydif = (it1.second.yy - it2->second.yy); double zdif = (it1.second.zz - it2->second.zz); - double vdif = std::max(std::abs(xdif), std::abs(ydif)); - vdif = std::max(vdif, std::abs(zdif)); - if (vdif > difmax1) { - difmax1 = vdif; - nameMax = it1.first; - } + double vdif = std::max(std::abs(xdif), std::abs(ydif)); + vdif = std::max(vdif, std::abs(zdif)); + if (vdif > difmax1) { + difmax1 = vdif; + nameMax = it1.first; + } if ((std::abs(xdif) > tol4) || (std::abs(ydif) > tol4) || (std::abs(zdif) > tol4)) { ++kount2; std::cout << it1.first << " x " << it1.second.xx << ":" << it2->second.xx << ":" << xdif << " y " diff --git a/Validation/Geometry/macros/MatBudgetVolume.C b/Validation/Geometry/macros/MatBudgetVolume.C index a793007f20012..aa51f0ae0f8aa 100644 --- a/Validation/Geometry/macros/MatBudgetVolume.C +++ b/Validation/Geometry/macros/MatBudgetVolume.C @@ -18,7 +18,7 @@ // Compares material budget plots from 2 different files // // etaPhiPlotComp4(filePreFix, tag, plot, ifEta, debug) -// Compares material budget plots from 4 different files: +// Compares material budget plots from 4 different files: // dddXML, dd4hepXML, dddDB, dd4hepDB // // filePreFix (std::string) Prefix to all 4 file names which will be followed @@ -66,13 +66,13 @@ void etaPhiPlotComp(TString fileName1 = "matbdg_run3.root", std::string plot = "intl", bool ifEta = true, std::string tag = "Run3", - std::string txt = "{DDD}/{DD4Hep}", + std::string txt = "{DDD}/{DD4Hep}", bool debug = false); -void etaPhiPlotComp4(std::string filePreFix = "files/matbdgRun3", - std::string tag = "pre6", - std::string plot = "radl", - bool ifEta = true, - bool debug = false); +void etaPhiPlotComp4(std::string filePreFix = "files/matbdgRun3", + std::string tag = "pre6", + std::string plot = "radl", + bool ifEta = true, + bool debug = false); void setStyle(); const int nlay = 13; @@ -242,7 +242,8 @@ void etaPhi2DPlot(TString fileName, std::string plot, bool drawLeg, double maxEt cc1->Modified(); } -void etaPhiPlotComp(TString fileName1, TString fileName2, std::string plot, bool ifEta, std::string tag, std::string txt, bool debug) { +void etaPhiPlotComp( + TString fileName1, TString fileName2, std::string plot, bool ifEta, std::string tag, std::string txt, bool debug) { setStyle(); gStyle->SetOptTitle(0); TFile *file1 = new TFile(fileName1); @@ -260,12 +261,12 @@ void etaPhiPlotComp(TString fileName1, TString fileName2, std::string plot, bool std::string ztit = "Eta"; char ytit[40]; if (plot == "radl") { - sprintf (ytit, "#frac%s for MB (X_{0})", txt.c_str()); + sprintf(ytit, "#frac%s for MB (X_{0})", txt.c_str()); } else if (plot == "step") { - sprintf (ytit, "#frac%s for MB (Step Length)", txt.c_str()); + sprintf(ytit, "#frac%s for MB (Step Length)", txt.c_str()); } else { plot = "intl"; - sprintf (ytit, "#frac%s for MB (#lambda)", txt.c_str()); + sprintf(ytit, "#frac%s for MB (#lambda)", txt.c_str()); } if (!ifEta) { xtit = "#phi"; @@ -313,19 +314,20 @@ void etaPhiPlotComp(TString fileName1, TString fileName2, std::string plot, bool } double temp1 = (rat > 1.0) ? 1.0 / rat : rat; double temp2 = (rat > 1.0) ? drt / (rat * rat) : drt; - double temp0 = (fabs(1.0 - temp1) / (temp2 * temp2)); + double temp0 = (fabs(1.0 - temp1) / (temp2 * temp2)); sumNum += temp0; sumDen += (1.0 / (temp2 * temp2)); - if (temp0 >= maxtmp) { - maxtmp = temp0; - maxDev = fabs(1.0 - temp1); - dmaxDev = temp2; - } + if (temp0 >= maxtmp) { + maxtmp = temp0; + maxDev = fabs(1.0 - temp1); + dmaxDev = temp2; + } } } sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; - std::cout << "Mean deviation for " << title[ii] << " " << sumNum << " +- " << sumDen << " Max " << maxDev << " +- " << dmaxDev << std::endl; + std::cout << "Mean deviation for " << title[ii] << " " << sumNum << " +- " << sumDen << " Max " << maxDev + << " +- " << dmaxDev << std::endl; if (xx.size() > 0) { TGraphErrors *graph = new TGraphErrors(xx.size(), &xx[0], &yy[0], &dx[0], &dy[0]); graph->SetLineColor(colorLay[ii]); @@ -421,71 +423,74 @@ void etaPhiPlotComp4(std::string filePreFix, std::string tag, std::string plot, int ii = nflayer[i] + j; sprintf(hname, "%s%s%s", plot.c_str(), ztit.c_str(), names[ii].c_str()); TProfile *prof[files]; - bool okf(true); - for (int k1 = 0; k1 < files; ++k1) { - dir[k1]->GetObject(hname, prof[k1]); - if (dir[k1] == nullptr) - okf = false; - } + bool okf(true); + for (int k1 = 0; k1 < files; ++k1) { + dir[k1]->GetObject(hname, prof[k1]); + if (dir[k1] == nullptr) + okf = false; + } if (okf) { int nb = prof[0]->GetNbinsX(); for (int k = 1; k <= nb; ++k) { xx0.push_back(prof[0]->GetBinLowEdge(k) + prof[0]->GetBinWidth(k)); - for (int k1 = 0; k1 < files; ++k1) { - yy0[k1].push_back(prof[k1]->GetBinContent(k)); - dy0[k1].push_back(prof[k1]->GetBinError(k)); - } + for (int k1 = 0; k1 < files; ++k1) { + yy0[k1].push_back(prof[k1]->GetBinContent(k)); + dy0[k1].push_back(prof[k1]->GetBinError(k)); + } } } } int ii = nflayer[i]; for (int k1 = 1; k1 < files; ++k1) { - std::vector xx, yy, dx, dy; - double sumNum(0), sumDen(0), maxtmp(0), maxDev(0), dmaxDev(0); - for (unsigned int k = 0; k < xx0.size(); ++k) { - if ((yy0[0][k] > 0) && (yy0[k1][k] > 0)) { - double rat = yy0[k1][k] / yy0[0][k]; - double drt = rat * sqrt((dy0[k1][k] / yy0[k1][k]) * (dy0[k1][k] / yy0[k1][k]) + (dy0[0][k] / yy0[0][k]) * (dy0[0][k] / yy0[0][k])); - xx.push_back(xx0[k]); - dx.push_back(0); - yy.push_back(rat); - dy.push_back(drt); - if (debug) { - std::cout << nametype[k1] << ":" << title[ii] << " [" << (xx.size() - 1) << "] " << xx0[k] << " Ratio " << rat << " +- " << drt << std::endl; + std::vector xx, yy, dx, dy; + double sumNum(0), sumDen(0), maxtmp(0), maxDev(0), dmaxDev(0); + for (unsigned int k = 0; k < xx0.size(); ++k) { + if ((yy0[0][k] > 0) && (yy0[k1][k] > 0)) { + double rat = yy0[k1][k] / yy0[0][k]; + double drt = rat * sqrt((dy0[k1][k] / yy0[k1][k]) * (dy0[k1][k] / yy0[k1][k]) + + (dy0[0][k] / yy0[0][k]) * (dy0[0][k] / yy0[0][k])); + xx.push_back(xx0[k]); + dx.push_back(0); + yy.push_back(rat); + dy.push_back(drt); + if (debug) { + std::cout << nametype[k1] << ":" << title[ii] << " [" << (xx.size() - 1) << "] " << xx0[k] << " Ratio " + << rat << " +- " << drt << std::endl; + } + double temp1 = (rat > 1.0) ? 1.0 / rat : rat; + double temp2 = (rat > 1.0) ? drt / (rat * rat) : drt; + double temp0 = (fabs(1.0 - temp1) / (temp2 * temp2)); + sumNum += temp0; + sumDen += (1.0 / (temp2 * temp2)); + if (temp0 >= maxtmp) { + maxtmp = temp0; + maxDev = fabs(1.0 - temp1); + dmaxDev = temp2; + } + } + } + sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; + sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; + std::cout << title[ii] << " in " << nametype[k1] << " Mean " << sumNum << " +- " << sumDen << " Max " << maxDev + << " +- " << dmaxDev << std::endl; + if (xx.size() > 0) { + TGraphErrors *graph = new TGraphErrors(xx.size(), &xx[0], &yy[0], &dx[0], &dy[0]); + graph->SetLineColor(colortype[k1]); + graph->SetFillColor(colorLay[ii]); + graph->SetMarkerStyle(styleLay[ii]); + if (k1 == 1) { + sprintf(titlex, "%s", title[ii].c_str()); + leg->AddEntry(graph, titlex, "lep"); + } + graphs.push_back(graph); + if (nb == 0) { + sprintf(hname, "%s%s%s", plot.c_str(), ztit.c_str(), names[0].c_str()); + TProfile *prof; + dir[0]->GetObject(hname, prof); + nb = prof->GetNbinsX(); + xlow = prof->GetBinLowEdge(1); + xhigh = prof->GetBinLowEdge(nb) + prof->GetBinWidth(nb); } - double temp1 = (rat > 1.0) ? 1.0 / rat : rat; - double temp2 = (rat > 1.0) ? drt / (rat * rat) : drt; - double temp0 = (fabs(1.0 - temp1) / (temp2 * temp2)); - sumNum += temp0; - sumDen += (1.0 / (temp2 * temp2)); - if (temp0 >= maxtmp) { - maxtmp = temp0; - maxDev = fabs(1.0 - temp1); - dmaxDev = temp2; - } - } - } - sumNum = (sumDen > 0) ? (sumNum / sumDen) : 0; - sumDen = (sumDen > 0) ? 1.0 / sqrt(sumDen) : 0; - std::cout << title[ii] << " in " << nametype[k1] << " Mean " << sumNum << " +- " << sumDen << " Max " << maxDev << " +- " << dmaxDev << std::endl; - if (xx.size() > 0) { - TGraphErrors *graph = new TGraphErrors(xx.size(), &xx[0], &yy[0], &dx[0], &dy[0]); - graph->SetLineColor(colortype[k1]); - graph->SetFillColor(colorLay[ii]); - graph->SetMarkerStyle(styleLay[ii]); - if (k1 == 1) { - sprintf(titlex, "%s", title[ii].c_str()); - leg->AddEntry(graph, titlex, "lep"); - } - graphs.push_back(graph); - if (nb == 0) { - sprintf(hname, "%s%s%s", plot.c_str(), ztit.c_str(), names[0].c_str()); - TProfile *prof; - dir[0]->GetObject(hname, prof); - nb = prof->GetNbinsX(); - xlow = prof->GetBinLowEdge(1); - xhigh = prof->GetBinLowEdge(nb) + prof->GetBinWidth(nb); - } } } } @@ -511,13 +516,13 @@ void etaPhiPlotComp4(std::string filePreFix, std::string tag, std::string plot, cc1->Modified(); double ymx = 0.68; for (int k1 = 1; k1 < files; ++k1) { - TPaveText* txt1 = new TPaveText(0.84, ymx - 0.03, 0.99, ymx, "blNDC"); - txt1->SetFillColor(0); - sprintf(fname, "%s", nametype[k1].c_str()); - txt1->AddText(fname); - ((TText*)txt1->GetListOfLines()->Last())->SetTextColor(colortype[k1]); - txt1->Draw(); - ymx -= 0.03; + TPaveText *txt1 = new TPaveText(0.84, ymx - 0.03, 0.99, ymx, "blNDC"); + txt1->SetFillColor(0); + sprintf(fname, "%s", nametype[k1].c_str()); + txt1->AddText(fname); + ((TText *)txt1->GetListOfLines()->Last())->SetTextColor(colortype[k1]); + txt1->Draw(); + ymx -= 0.03; } } }