Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Run3-sim107 Modify the test for SimHit comparisons #36763

Merged
merged 5 commits into from Jan 23, 2022
Merged

Run3-sim107 Modify the test for SimHit comparisons #36763

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
f5e794c
Modify the test for SimHit comparisons
Jan 20, 2022
2690ad5
Code check + modified cfg
Jan 20, 2022
c2f2fc8
Small changes
Jan 21, 2022
02075c6
Fix the x-range
Jan 22, 2022
492d98b
Code check
Jan 22, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
132 changes: 132 additions & 0 deletions SimG4CMS/Calo/macros/MakeHitStudyPlots.C
View file
Open in desktop
@@ -0,0 +1,132 @@
#include <TCanvas.h>
#include <TChain.h>
#include <TFile.h>
#include <TFitResult.h>
#include <TFitResultPtr.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TLegend.h>
#include <TPaveStats.h>
#include <TPaveText.h>
#include <TProfile.h>
#include <TROOT.h>
#include <TStyle.h>
#include <vector>
#include <string>
#include <iomanip>
#include <iostream>
#include <fstream>

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) {
std::string names[18] = {"Edep",
"EdepEM",
"EdepHad",
"EdepTk",
"Etot",
"EtotG",
"Hit",
"HitHigh",
"HitLow",
"HitMu",
"HitTk",
"Time",
"TimeAll",
"TimeTk",
"EneInc",
"EtaInc",
"PhiInc",
"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};

gStyle->SetCanvasBorderMode(0);
gStyle->SetCanvasColor(kWhite);
gStyle->SetPadColor(kWhite);
gStyle->SetFillColor(kWhite);
gStyle->SetOptStat(1110);
TFile* file[2];
int nfile(0);
std::string tag(""), tags[2];
if (file1 != "") {
file[nfile] = new TFile(file1.c_str());
if (file[nfile]) {
tags[nfile] = tag1;
++nfile;
tag += tag1;
}
}
if (file2 != "") {
file[nfile] = new TFile(file2.c_str());
if (file[nfile]) {
tags[nfile] = tag2;
++nfile;
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());
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());
}
TH1D* hist = static_cast<TH1D*>(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);
}
}
}
}
}
78 changes: 58 additions & 20 deletions SimG4CMS/Calo/plugins/CaloSimHitStudy.cc
View file
Open in desktop
Expand Up @@ -33,6 +33,8 @@
#include <map>
#include <string>

//#define EDM_ML_DEBUG

class CaloSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
public:
CaloSimHitStudy(const edm::ParameterSet& ps);
Expand All @@ -53,17 +55,17 @@ class CaloSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::on
std::string g4Label, hitLab[4];
edm::EDGetTokenT<edm::HepMCProduct> tok_evt_;
edm::EDGetTokenT<edm::PCaloHitContainer> toks_calo_[4];
edm::EDGetTokenT<edm::PSimHitContainer> toks_track_[3];
edm::EDGetTokenT<edm::PSimHitContainer> toks_track_[4];
edm::EDGetTokenT<edm::PSimHitContainer> toks_tkHigh_[6];
edm::EDGetTokenT<edm::PSimHitContainer> toks_tkLow_[6];
std::string muonLab[3], tkHighLab[6], tkLowLab[6];
std::string muonLab[4], tkHighLab[6], tkLowLab[6];
double tmax_, eMIP_;
bool storeRL_, testNumber_;

TH1F *hit_[9], *time_[9], *edepEM_[9], *edepHad_[9], *edep_[9];
TH1F *etot_[9], *etotg_[9], *timeAll_[9], *hitMu, *hitHigh;
TH1F *hitLow, *eneInc_, *etaInc_, *phiInc_, *ptInc_;
TH1F *hitTk_[15], *edepTk_[15], *tofTk_[15], *edepC_[9], *edepT_[9];
TH1F *hitTk_[16], *edepTk_[16], *tofTk_[16], *edepC_[9], *edepT_[9];
};

CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
Expand All @@ -81,11 +83,12 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
tmax_ = ps.getUntrackedParameter<double>("TimeCut", 100.0);
eMIP_ = ps.getUntrackedParameter<double>("MIPCut", 0.70);
storeRL_ = ps.getUntrackedParameter<bool>("StoreRL", false);
testNumber_ = ps.getUntrackedParameter<bool>("TestNumbering", false);
testNumber_ = ps.getUntrackedParameter<bool>("TestNumbering", true);

muonLab[0] = "MuonRPCHits";
muonLab[1] = "MuonCSCHits";
muonLab[2] = "MuonDTHits";
muonLab[3] = "MuonGEMHits";
tkHighLab[0] = "TrackerHitsPixelBarrelHighTof";
tkHighLab[1] = "TrackerHitsPixelEndcapHighTof";
tkHighLab[2] = "TrackerHitsTECHighTof";
Expand All @@ -103,7 +106,7 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
for (unsigned i = 0; i != 4; i++)
toks_calo_[i] = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hitLab[i]));

for (unsigned i = 0; i != 3; i++)
for (unsigned i = 0; i != 4; i++)
toks_track_[i] = consumes<edm::PSimHitContainer>(edm::InputTag(g4Label, muonLab[i]));

for (unsigned i = 0; i != 6; i++) {
Expand Down Expand Up @@ -137,11 +140,15 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
phiInc_ = tfile->make<TH1F>("PhiInc", title, 200, -3.1415926, 3.1415926);
phiInc_->GetXaxis()->SetTitle(title);
phiInc_->GetYaxis()->SetTitle("Events");
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for incident particle";
#endif
std::string dets[9] = {"EB", "EB(APD)", "EB(ATJ)", "EE", "ES", "HB", "HE", "HO", "HF"};
double nhcMax[9] = {40000., 2000., 2000., 40000., 10000., 10000., 10000., 2000., 10000.};
for (int i = 0; i < 9; i++) {
sprintf(name, "Hit%d", i);
sprintf(title, "Number of hits in %s", dets[i].c_str());
hit_[i] = tfile->make<TH1F>(name, title, 1000, 0., 20000.);
hit_[i] = tfile->make<TH1F>(name, title, 1000, 0., nhcMax[i]);
hit_[i]->GetXaxis()->SetTitle(title);
hit_[i]->GetYaxis()->SetTitle("Events");
sprintf(name, "Time%d", i);
Expand Down Expand Up @@ -185,6 +192,9 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
etotg_[i]->GetXaxis()->SetTitle(title);
etotg_[i]->GetYaxis()->SetTitle("Events");
}
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for first level of Calorimeter";
#endif
std::string detx[9] = {"EB/EE (MIP)",
"HB/HE (MIP)",
"HB/HE/HO (MIP)",
Expand All @@ -209,16 +219,22 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
edepT_[i]->GetXaxis()->SetTitle(title);
edepT_[i]->GetYaxis()->SetTitle("Events");
}
hitLow = tfile->make<TH1F>("HitLow", "Number of hits in Track (Low)", 1000, 0, 10000.);
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for second level of Calorimeter";
#endif
hitLow = tfile->make<TH1F>("HitLow", "Number of hits in Track (Low)", 1000, 0, 20000.);
hitLow->GetXaxis()->SetTitle("Number of hits in Track (Low)");
hitLow->GetYaxis()->SetTitle("Events");
hitHigh = tfile->make<TH1F>("HitHigh", "Number of hits in Track (High)", 1000, 0, 10000.);
hitHigh = tfile->make<TH1F>("HitHigh", "Number of hits in Track (High)", 1000, 0, 5000.);
hitHigh->GetXaxis()->SetTitle("Number of hits in Track (High)");
hitHigh->GetYaxis()->SetTitle("Events");
hitMu = tfile->make<TH1F>("HitMu", "Number of hits in Track (Muon)", 1000, 0, 5000.);
hitMu = tfile->make<TH1F>("HitMu", "Number of hits in Track (Muon)", 1000, 0, 2000.);
hitMu->GetXaxis()->SetTitle("Number of hits in Muon");
hitMu->GetYaxis()->SetTitle("Events");
std::string dett[15] = {"Pixel Barrel (High)",
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for general tracking hits";
#endif
std::string dett[16] = {"Pixel Barrel (High)",
"Pixel Endcap (High)",
"TEC (High)",
"TIB (High)",
Expand All @@ -232,11 +248,14 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
"TOB (Low)",
"RPC",
"CSC",
"DT"};
for (int i = 0; i < 15; i++) {
"DT",
"GEM"};
double nhtMax[16] = {
500., 500., 1000., 1000., 500., 1000., 5000., 2000., 10000., 5000., 2000., 5000., 500., 1000., 1000., 500.};
for (int i = 0; i < 16; i++) {
sprintf(name, "HitTk%d", i);
sprintf(title, "Number of hits in %s", dett[i].c_str());
hitTk_[i] = tfile->make<TH1F>(name, title, 1000, 0., 1000.);
hitTk_[i] = tfile->make<TH1F>(name, title, 1000, 0., nhtMax[i]);
hitTk_[i]->GetXaxis()->SetTitle(title);
hitTk_[i]->GetYaxis()->SetTitle("Events");
sprintf(name, "TimeTk%d", i);
Expand All @@ -246,10 +265,14 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) {
tofTk_[i]->GetYaxis()->SetTitle("Hits");
sprintf(name, "EdepTk%d", i);
sprintf(title, "Energy deposit (GeV) in %s", dett[i].c_str());
edepTk_[i] = tfile->make<TH1F>(name, title, 5000, 0., 10.);
double ymax = (i < 12) ? 0.25 : 0.005;
edepTk_[i] = tfile->make<TH1F>(name, title, 250, 0., ymax);
edepTk_[i]->GetXaxis()->SetTitle(title);
edepTk_[i]->GetYaxis()->SetTitle("Hits");
}
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for SimHit objects";
#endif
}

void CaloSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
Expand All @@ -264,7 +287,7 @@ void CaloSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descripti
desc.addUntracked<double>("TimeCut", 100.0);
desc.addUntracked<double>("MIPCut", 0.70);
desc.addUntracked<bool>("StoreRL", false);
desc.addUntracked<bool>("TestNumbering", false);
desc.addUntracked<bool>("TestNumbering", true);
descriptions.add("CaloSimHitStudy", desc);
}

Expand Down Expand Up @@ -295,8 +318,9 @@ void CaloSimHitStudy::analyze(edm::Event const& e, edm::EventSetup const&) {
e.getByToken(toks_calo_[i], hitsCalo);
if (hitsCalo.isValid())
getHits = true;
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Input flags Hits " << getHits;

#endif
if (getHits) {
std::vector<PCaloHit> caloHits;
caloHits.insert(caloHits.end(), hitsCalo->begin(), hitsCalo->end());
Expand All @@ -306,15 +330,17 @@ void CaloSimHitStudy::analyze(edm::Event const& e, edm::EventSetup const&) {
eeHits.insert(eeHits.end(), hitsCalo->begin(), hitsCalo->end());
else if (i == 3)
hcHits.insert(hcHits.end(), hitsCalo->begin(), hitsCalo->end());
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Hit buffer " << caloHits.size();
#endif
analyzeHits(caloHits, i);
}
}
analyzeHits(ebHits, eeHits, hcHits);

std::vector<PSimHit> muonHits;
edm::Handle<edm::PSimHitContainer> hitsTrack;
for (int i = 0; i < 3; i++) {
for (int i = 0; i < 4; i++) {
e.getByToken(toks_track_[i], hitsTrack);
if (hitsTrack.isValid()) {
muonHits.insert(muonHits.end(), hitsTrack->begin(), hitsTrack->end());
Expand Down Expand Up @@ -347,7 +373,10 @@ void CaloSimHitStudy::analyze(edm::Event const& e, edm::EventSetup const&) {

void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
int nHit = hits.size();
int nHB = 0, nHE = 0, nHO = 0, nHF = 0, nEB = 0, nEBAPD = 0, nEBATJ = 0, nEE = 0, nES = 0, nBad = 0;
int nHB(0), nHE(0), nHO(0), nHF(0), nEB(0), nEBAPD(0), nEBATJ(0), nEE(0), nES(0);
#ifdef EDM_ML_DEBUG
int nBad(0);
#endif
std::map<unsigned int, double> hitMap;
std::vector<double> etot(9, 0), etotG(9, 0);
for (int i = 0; i < nHit; i++) {
Expand Down Expand Up @@ -390,8 +419,10 @@ void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
nEE++;
else if (idx == 4)
nES++;
#ifdef EDM_ML_DEBUG
else
nBad++;
#endif
if (indx >= 0 && indx < 3) {
etot[idx] += edep;
if (time < 100)
Expand Down Expand Up @@ -428,7 +459,9 @@ void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
if (time < 100)
etotG[idx] += edep;
} else {
#ifdef EDM_ML_DEBUG
nBad++;
#endif
}
}
}
Expand Down Expand Up @@ -457,9 +490,11 @@ void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
}
}

#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::analyzeHits: EB " << nEB << ", " << nEBAPD << ", " << nEBATJ
<< " EE " << nEE << " ES " << nES << " HB " << nHB << " HE " << nHE << " HO " << nHO
<< " HF " << nHF << " Bad " << nBad << " All " << nHit << " Reduced " << hitMap.size();
#endif
std::map<unsigned int, double>::const_iterator it = hitMap.begin();
for (; it != hitMap.end(); it++) {
double time = it->second;
Expand Down Expand Up @@ -511,16 +546,18 @@ void CaloSimHitStudy::analyzeHits(edm::Handle<edm::PSimHitContainer>& hits, int
label = tkHighLab[indx];
else if (indx >= 6 && indx < 12)
label = tkLowLab[indx - 6];
else if (indx >= 12 && indx < 15)
else if (indx >= 12 && indx < 16)
label = muonLab[indx - 12];
for (ihit = hits->begin(); ihit != hits->end(); ihit++) {
edepTk_[indx]->Fill(ihit->energyLoss());
tofTk_[indx]->Fill(ihit->timeOfFlight());
nHit++;
}
hitTk_[indx]->Fill(float(nHit));
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::analyzeHits: for " << label << " Index " << indx << " # of Hits "
<< nHit;
#endif
}

void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& ebHits,
Expand Down Expand Up @@ -570,12 +607,13 @@ void CaloSimHitStudy::analyzeHits(std::vector<PCaloHit>& ebHits,
double edepET = edepEBT + edepEET;
double edepHC = edepH + edepHO;
double edepHCT = edepHT + edepHOT;
#ifdef EDM_ML_DEBUG
edm::LogVerbatim("HitStudy") << "CaloSimHitStudy::energy in EB " << edepEB << " (" << edepEBT << ") from "
<< ebHits.size() << " hits; "
<< " energy in EE " << edepEE << " (" << edepEET << ") from " << eeHits.size()
<< " hits; energy in HC " << edepH << ", " << edepHO << " (" << edepHT << ", " << edepHOT
<< ") from " << hcHits.size() << " hits";

#endif
edepC_[6]->Fill(edepE);
edepT_[6]->Fill(edepET);
edepC_[7]->Fill(edepH);
Expand Down