From 9202783f6cb81eaf725cc93e332961b67c6c20f3 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sun, 6 Nov 2022 09:51:28 +0100 Subject: [PATCH 1/4] Modify the analysis code of SimHit Results --- SimG4CMS/Calo/plugins/CaloSimHitStudy.cc | 131 ++++++++++++++-------- SimG4CMS/Calo/python/calostep_cfi.py | 53 +++++++++ SimG4CMS/Calo/test/python/minbias.py | 135 +++++++++++++++++++++++ SimG4CMS/Calo/test/python/ttbar.py | 131 ++++++++++++++++++++++ 4 files changed, 401 insertions(+), 49 deletions(-) create mode 100644 SimG4CMS/Calo/python/calostep_cfi.py create mode 100644 SimG4CMS/Calo/test/python/minbias.py create mode 100644 SimG4CMS/Calo/test/python/ttbar.py diff --git a/SimG4CMS/Calo/plugins/CaloSimHitStudy.cc b/SimG4CMS/Calo/plugins/CaloSimHitStudy.cc index 67e8a67a1259a..6c7c8a7410b17 100644 --- a/SimG4CMS/Calo/plugins/CaloSimHitStudy.cc +++ b/SimG4CMS/Calo/plugins/CaloSimHitStudy.cc @@ -15,7 +15,10 @@ #include "FWCore/ServiceRegistry/interface/Service.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" #include "DataFormats/HcalDetId/interface/HcalDetId.h" +#include "DataFormats/HcalDetId/interface/HcalSubdetector.h" #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h" #include "SimDataFormats/CaloHit/interface/PCaloHit.h" @@ -25,6 +28,11 @@ #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h" #include "SimG4CMS/Calo/interface/CaloHitID.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" + #include #include @@ -57,12 +65,16 @@ class CaloSimHitStudy : public edm::one::EDAnalyzer hitLab_; const double maxEnergy_, tmax_, eMIP_; const bool storeRL_, testNumber_; - edm::EDGetTokenT tok_evt_; - std::vector> toks_calo_; + const edm::ESGetToken tokGeom_; + const edm::EDGetTokenT tok_evt_; + const std::vector> toks_calo_; std::vector> toks_track_; std::vector> toks_tkHigh_; std::vector> toks_tkLow_; + const CaloGeometry* caloGeometry_; + const HcalGeometry* hcalGeom_; + const std::vector muonLab_ = {"MuonRPCHits", "MuonCSCHits", "MuonDTHits", "MuonGEMHits"}; const std::vector tkHighLab_ = {"TrackerHitsPixelBarrelHighTof", "TrackerHitsPixelEndcapHighTof", @@ -77,10 +89,13 @@ class CaloSimHitStudy : public edm::one::EDAnalyzer("TimeCut", 100.0)), eMIP_(ps.getUntrackedParameter("MIPCut", 0.70)), storeRL_(ps.getUntrackedParameter("StoreRL", false)), - testNumber_(ps.getUntrackedParameter("TestNumbering", true)) { + testNumber_(ps.getUntrackedParameter("TestNumbering", true)), + tokGeom_(esConsumes()), + tok_evt_(consumes(edm::InputTag(ps.getUntrackedParameter("SourceLabel", "VtxSmeared")))), + toks_calo_{edm::vector_transform(hitLab_, + [this](const std::string& name) { + return consumes(edm::InputTag{g4Label_, name}); + })} { usesResource(TFileService::kSharedResource); - tok_evt_ = - consumes(edm::InputTag(ps.getUntrackedParameter("SourceLabel", "VtxSmeared"))); - for (unsigned i = 0; i != hitLab_.size(); i++) - toks_calo_.emplace_back(consumes(edm::InputTag(g4Label_, hitLab_[i]))); for (unsigned i = 0; i != muonLab_.size(); i++) toks_track_.emplace_back(consumes(edm::InputTag(g4Label_, muonLab_[i]))); for (unsigned i = 0; i != tkHighLab_.size(); i++) @@ -133,9 +150,9 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) #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++) { + std::string dets[nCalo_] = {"EB", "EB(APD)", "EB(ATJ)", "EE", "ES", "HB", "HE", "HO", "HF"}; + double nhcMax[nCalo_] = {40000., 2000., 2000., 40000., 10000., 10000., 10000., 2000., 10000.}; + for (int i = 0; i < nCalo_; i++) { sprintf(name, "Hit%d", i); sprintf(title, "Number of hits in %s", dets[i].c_str()); hit_[i] = tfile->make(name, title, 1000, 0., nhcMax[i]); @@ -181,6 +198,16 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) etotg_[i] = tfile->make(name, title, 5000, 0., ymax); etotg_[i]->GetXaxis()->SetTitle(title); etotg_[i]->GetYaxis()->SetTitle("Events"); + sprintf(name, "eta%d", i); + sprintf(title, "#eta of hit point in %s", dets[i].c_str()); + eta_[i] = tfile->make(name, title, 100, -5.0, 5.0); + eta_[i]->GetXaxis()->SetTitle(title); + eta_[i]->GetYaxis()->SetTitle("Hits"); + sprintf(name, "phi%d", i); + sprintf(title, "#phi of hit point in %s", dets[i].c_str()); + phi_[i] = tfile->make(name, title, 100, -M_PI, M_PI); + phi_[i]->GetXaxis()->SetTitle(title); + phi_[i]->GetYaxis()->SetTitle("Hits"); } #ifdef EDM_ML_DEBUG edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for first level of Calorimeter"; @@ -224,25 +251,25 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) #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)", - "TID (High)", - "TOB (High)", - "Pixel Barrel (Low)", - "Pixel Endcap (Low)", - "TEC (Low)", - "TIB (Low)", - "TID (Low)", - "TOB (Low)", - "RPC", - "CSC", - "DT", - "GEM"}; - double nhtMax[16] = { + std::string dett[nTrack_] = {"Pixel Barrel (High)", + "Pixel Endcap (High)", + "TEC (High)", + "TIB (High)", + "TID (High)", + "TOB (High)", + "Pixel Barrel (Low)", + "Pixel Endcap (Low)", + "TEC (Low)", + "TIB (Low)", + "TID (Low)", + "TOB (Low)", + "RPC", + "CSC", + "DT", + "GEM"}; + double nhtMax[nTrack_] = { 500., 500., 1000., 1000., 500., 1000., 5000., 2000., 10000., 5000., 2000., 5000., 500., 1000., 1000., 500.}; - for (int i = 0; i < 16; i++) { + for (int i = 0; i < nTrack_; i++) { sprintf(name, "HitTk%d", i); sprintf(title, "Number of hits in %s", dett[i].c_str()); hitTk_[i] = tfile->make(name, title, 1000, 0., nhtMax[i]); @@ -279,9 +306,12 @@ void CaloSimHitStudy::fillDescriptions(edm::ConfigurationDescriptions& descripti descriptions.add("CaloSimHitStudy", desc); } -void CaloSimHitStudy::analyze(edm::Event const& e, edm::EventSetup const&) { +void CaloSimHitStudy::analyze(edm::Event const& e, edm::EventSetup const& set) { edm::LogVerbatim("HitStudy") << "CaloSimHitStudy:Run = " << e.id().run() << " Event = " << e.id().event(); + caloGeometry_ = &set.getData(tokGeom_); + hcalGeom_ = static_cast(caloGeometry_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel)); + const edm::Handle EvtHandle = e.getHandle(tok_evt_); const HepMC::GenEvent* myGenEvent = EvtHandle->GetEvent(); @@ -363,23 +393,23 @@ void CaloSimHitStudy::analyzeHits(std::vector& hits, int indx) { int nBad(0); #endif std::map hitMap; - std::vector etot(9, 0), etotG(9, 0); + std::vector etot(nCalo_, 0), etotG(nCalo_, 0); for (int i = 0; i < nHit; i++) { double edep = hits[i].energy(); double time = hits[i].time(); - unsigned int id_ = hits[i].id(); + unsigned int id = hits[i].id(); double edepEM = hits[i].energyEM(); double edepHad = hits[i].energyHad(); if (indx == 0) { int dep = (hits[i].depth()) & PCaloHit::kEcalDepthIdMask; if (dep == 1) - id_ |= 0x20000; + id |= 0x20000; else if (dep == 2) - id_ |= 0x40000; + id |= 0x40000; } - std::map::const_iterator it = hitMap.find(id_); + std::map::const_iterator it = hitMap.find(id); if (it == hitMap.end()) { - hitMap.insert(std::pair(id_, time)); + hitMap.insert(std::pair(id, time)); } int idx = -1; if (indx != 3) { @@ -417,10 +447,9 @@ void CaloSimHitStudy::analyzeHits(std::vector& hits, int indx) { int subdet(0); if (testNumber_) { int ieta(0), phi(0), z(0), lay(0), depth(0); - HcalTestNumbering::unpackHcalIndex(id_, subdet, z, depth, ieta, phi, lay); + HcalTestNumbering::unpackHcalIndex(id, subdet, z, depth, ieta, phi, lay); } else { - HcalDetId id = HcalDetId(id_); - subdet = id.subdet(); + subdet = HcalDetId(id).subdet(); } if (subdet == static_cast(HcalBarrel)) { idx = indx + 2; @@ -483,13 +512,16 @@ void CaloSimHitStudy::analyzeHits(std::vector& hits, int indx) { std::map::const_iterator it = hitMap.begin(); for (; it != hitMap.end(); it++) { double time = it->second; - unsigned int id_ = (it->first); + GlobalPoint point; + DetId id(it->first); + if (indx != 2) + point = caloGeometry_->getPosition(id); int idx = -1; if (indx < 3) { if (indx == 0) { - if ((id_ & 0x20000) != 0) + if ((id & 0x20000) != 0) idx = indx + 1; - else if ((id_ & 0x40000) != 0) + else if ((id & 0x40000) != 0) idx = indx + 1; else idx = indx; @@ -502,10 +534,9 @@ void CaloSimHitStudy::analyzeHits(std::vector& hits, int indx) { int idx(-1), subdet(0); if (testNumber_) { int ieta(0), phi(0), z(0), lay(0), depth(0); - HcalTestNumbering::unpackHcalIndex(id_, subdet, z, depth, ieta, phi, lay); + HcalTestNumbering::unpackHcalIndex(id.rawId(), subdet, z, depth, ieta, phi, lay); } else { - HcalDetId id = HcalDetId(id_); - subdet = id.subdet(); + subdet = HcalDetId(id).subdet(); } if (subdet == static_cast(HcalBarrel)) { idx = indx + 2; @@ -518,6 +549,8 @@ void CaloSimHitStudy::analyzeHits(std::vector& hits, int indx) { } if (idx > 0) { timeAll_[idx]->Fill(time); + eta_[idx]->Fill(point.eta()); + phi_[idx]->Fill(point.phi()); } } } @@ -531,7 +564,7 @@ void CaloSimHitStudy::analyzeHits(const edm::Handle& hits label = tkHighLab_[indx]; else if (indx >= 6 && indx < 12) label = tkLowLab_[indx - 6]; - else if (indx >= 12 && indx < 16) + else if (indx >= 12 && indx < nTrack_) label = muonLab_[indx - 12]; for (ihit = hits->begin(); ihit != hits->end(); ihit++) { edepTk_[indx]->Fill(ihit->energyLoss()); diff --git a/SimG4CMS/Calo/python/calostep_cfi.py b/SimG4CMS/Calo/python/calostep_cfi.py new file mode 100644 index 0000000000000..6371266973dcc --- /dev/null +++ b/SimG4CMS/Calo/python/calostep_cfi.py @@ -0,0 +1,53 @@ +import FWCore.ParameterSet.Config as cms +from SimG4Core.Configuration.SimG4Core_cff import * + +g4SimHits.Watchers = cms.VPSet(cms.PSet( + type = cms.string('CaloSteppingAction'), + CaloSteppingAction = cms.PSet( + EBSDNames = cms.vstring('EBRY'), + EESDNames = cms.vstring('EFRY'), + HCSDNames = cms.vstring('HBS','HES','HTS'), + AllSteps = cms.int32(100), + SlopeLightYield = cms.double(0.02), + BirkC1EC = cms.double(0.03333), + BirkSlopeEC = cms.double(0.253694), + BirkCutEC = cms.double(0.1), + BirkC1HC = cms.double(0.0052), + BirkC2HC = cms.double(0.142), + BirkC3HC = cms.double(1.75), + HitCollNames = cms.vstring('EcalHitsEB1','EcalHitsEE1', + 'HcalHits1'), + EtaTable = cms.vdouble(0.000, 0.087, 0.174, 0.261, 0.348, + 0.435, 0.522, 0.609, 0.696, 0.783, + 0.870, 0.957, 1.044, 1.131, 1.218, + 1.305, 1.392, 1.479, 1.566, 1.653, + 1.740, 1.830, 1.930, 2.043, 2.172, + 2.322, 2.500, 2.650, 2.868, 3.000), + PhiBin = cms.vdouble(5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, + 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, + 5.0, 5.0, 5.0, 5.0, 5.0, 5.0,10.0, + 10.0,10.0,10.0,10.0,10.0,10.0,10.0, + 10.0), + PhiOffset = cms.vdouble( 0.0, 0.0, 0.0,10.0), + EtaMin = cms.vint32(1, 16, 29, 1), + EtaMax = cms.vint32(16, 29, 41, 15), + EtaHBHE = cms.int32(16), + DepthHBHE = cms.vint32(2,4), + Depth29Max = cms.int32(3), + RMinHO = cms.double(3800), + ZHO = cms.vdouble(0,1255,1418,3928,4100,6610), + Eta1 = cms.untracked.vint32(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 1, 1, 1, 1, 4, 4), + Eta15 = cms.untracked.vint32(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, + 1, 1, 1, 2, 2, 2, 2, 4, 4), + Eta16 = cms.untracked.vint32(1, 1, 2, 2, 2, 2, 2, 2, 2, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4), + Eta17 = cms.untracked.vint32(2, 2, 2, 2, 2, 2, 2, 2, 3, 3, + 3, 3, 3, 3, 3, 3, 3, 3, 3), + Eta18 = cms.untracked.vint32(1, 2, 2, 2, 3, 3, 3, 3, 4, 4, + 4, 5, 5, 5, 5, 5, 5, 5, 5), + Eta19 = cms.untracked.vint32(1, 1, 2, 2, 2, 3, 3, 3, 4, 4, + 4, 5, 5, 5, 5, 6, 6, 6, 6), + Eta26 = cms.untracked.vint32(1, 1, 2, 2, 3, 3, 4, 4, 5, 5, + 5, 6, 6, 6, 6, 7, 7, 7, 7), + ))) diff --git a/SimG4CMS/Calo/test/python/minbias.py b/SimG4CMS/Calo/test/python/minbias.py new file mode 100644 index 0000000000000..897d0226ea4e7 --- /dev/null +++ b/SimG4CMS/Calo/test/python/minbias.py @@ -0,0 +1,135 @@ +############################################################################### +# Way to use this: +# cmsRun minbias.py geometry=2016 +# Options for geometry 2016, 2017, 2018, 2021, legacy +# +############################################################################### +import FWCore.ParameterSet.Config as cms +import os, sys, imp, re +import FWCore.ParameterSet.VarParsing as VarParsing + +#################################################################### +### SETUP OPTIONS +options = VarParsing.VarParsing('standard') +options.register('geometry', + "2021", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "geometry of operations: 2016, 2017, 2018, 2021") + +### get and parse the command line arguments + +options.parseArguments() + +print(options) + +#################################################################### +# Use the options +histFile = "minbias" + options.geometry + ".root" +outFile = "simevent_minbias" + options.geometry + ".root" + +if (options.geometry == "2016"): + from Configuration.Eras.Era_Run2_2016_cff import Run2_2016 + process = cms.Process('Sim',Run2_2016) + geomFile = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff" + globalTag = "auto:run2_mc" +elif (options.geometry == "2017"): + from Configuration.Eras.Era_Run2_2017_cff import Run2_2017 + process = cms.Process('Sim',Run2_2017) + geomFile = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff" + globalTag = "auto:phase1_2017_realistic" +elif (options.geometry == "2018"): + from Configuration.Eras.Era_Run2_2018_cff import Run2_2018 + process = cms.Process('Sim',Run2_2018) + geomFile = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff" + globalTag = "auto:phase1_2018_realistic" +elif (options.geometry == "2021"): + from Configuration.Eras.Era_Run3_DDD_cff import Run3_DDD + process = cms.Process('Sim',Run3_DDD) + geomFile = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff" + globalTag = "auto:phase1_2022_realistic" +else: + process = cms.Process('Sim') + geomFile = "Configuration.Geometry.GeometryExtendedReco_cff" + globalTag = "auto:run1_mc" + +print("Geometry file: ", geomFile) +print("Hist file: ", histFile) +print("Output file: ", outFile) +print("Gobal Tag: ", globalTag) + +process.load("SimG4CMS.Calo.PythiaMinBias_cfi") +process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load("IOMC.EventVertexGenerators.VtxSmearedGauss_cfi") +process.load(geomFile) +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.EventContent.EventContent_cff") +process.load('Configuration.StandardSequences.Generator_cff') +process.load('Configuration.StandardSequences.SimIdeal_cff') +process.load("SimG4CMS.Calo.CaloSimHitStudy_cfi") +process.load("SimG4CMS.Calo.hcalSimHitDump_cfi") + +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, globalTag, '') + +if 'MessageLogger' in process.__dict__: + process.MessageLogger.G4cerr=dict() + process.MessageLogger.HitStudy=dict() + process.MessageLogger.HcalSim=dict() + +process.source = cms.Source("EmptySource") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(1000) +) + +process.Timing = cms.Service("Timing") + +process.SimpleMemoryCheck = cms.Service("SimpleMemoryCheck", + oncePerEventMode = cms.untracked.bool(True), + showMallocInfo = cms.untracked.bool(True), + dump = cms.untracked.bool(True), + ignoreTotal = cms.untracked.int32(1) +) + +process.load("IOMC.RandomEngine.IOMC_cff") +process.RandomNumberGeneratorService.generator.initialSeed = 456789 +process.RandomNumberGeneratorService.g4SimHits.initialSeed = 9876 +process.RandomNumberGeneratorService.VtxSmeared.initialSeed = 123456789 +process.rndmStore = cms.EDProducer("RandomEngineStateProducer") + +process.TFileService = cms.Service("TFileService", + fileName = cms.string(histFile) +) + +# Event output +process.output = cms.OutputModule("PoolOutputModule", + process.FEVTSIMEventContent, + fileName = cms.untracked.string(outFile) +) + +process.generation_step = cms.Path(process.pgen) +process.simulation_step = cms.Path(process.psim) +process.analysis_step = cms.Path(process.CaloSimHitStudy+process.hcalSimHitDump) +process.out_step = cms.EndPath(process.output) + +process.generator.pythiaHepMCVerbosity = False +process.generator.pythiaPylistVerbosity = 0 +process.g4SimHits.Physics.type = 'SimG4Core/Physics/FTFP_BERT_EMM' +process.g4SimHits.HCalSD.TestNumberingScheme = False +process.CaloSimHitStudy.TestNumbering = False +process.hcalSimHitDump.MaxEvent = 5 + +# Schedule definition +process.schedule = cms.Schedule(process.generation_step, + process.simulation_step, + process.analysis_step, +# process.out_step + ) + +# filter all path with the production filter sequence +for path in process.paths: + getattr(process,path)._seq = process.generator * getattr(process,path)._seq + diff --git a/SimG4CMS/Calo/test/python/ttbar.py b/SimG4CMS/Calo/test/python/ttbar.py new file mode 100644 index 0000000000000..e9911023635e4 --- /dev/null +++ b/SimG4CMS/Calo/test/python/ttbar.py @@ -0,0 +1,131 @@ +############################################################################### +# Way to use this: +# cmsRun ttbar.py geometry=2016 +# Options for geometry 2016, 2017, 2018, 2021, legacy +# +############################################################################### +import FWCore.ParameterSet.Config as cms +import os, sys, imp, re +import FWCore.ParameterSet.VarParsing as VarParsing + +#################################################################### +### SETUP OPTIONS +options = VarParsing.VarParsing('standard') +options.register('geometry', + "2021", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "geometry of operations: 2016, 2017, 2018, 2021") + +### get and parse the command line arguments + +options.parseArguments() + +print(options) + +#################################################################### +# Use the options +histFile = "ttbar" + options.geometry + ".root" +outFile = "simevent_ttbar" + options.geometry + ".root" + +if (options.geometry == "2016"): + from Configuration.Eras.Era_Run2_2016_cff import Run2_2016 + process = cms.Process('Sim',Run2_2016) + geomFile = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff" + globalTag = "auto:run2_mc" +elif (options.geometry == "2017"): + from Configuration.Eras.Era_Run2_2017_cff import Run2_2017 + process = cms.Process('Sim',Run2_2017) + geomFile = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff" + globalTag = "auto:phase1_2017_realistic" +elif (options.geometry == "2018"): + from Configuration.Eras.Era_Run2_2018_cff import Run2_2018 + process = cms.Process('Sim',Run2_2018) + geomFile = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff" + globalTag = "auto:phase1_2018_realistic" +elif (options.geometry == "2021"): + from Configuration.Eras.Era_Run3_DDD_cff import Run3_DDD + process = cms.Process('Sim',Run3_DDD) + geomFile = "Configuration.Geometry.GeometryExtended" + options.geometry + "Reco_cff" + globalTag = "auto:phase1_2022_realistic" +else: + process = cms.Process('Sim') + geomFile = "Configuration.Geometry.GeometryExtendedReco_cff" + globalTag = "auto:run1_mc" + +print("Geometry file: ", geomFile) +print("Hist file: ", histFile) +print("Output file: ", outFile) +print("Gobal Tag: ", globalTag) + +process.load("SimG4CMS.Calo.PythiaTT_cfi") +process.load("SimGeneral.HepPDTESSource.pythiapdt_cfi") +process.load('FWCore.MessageService.MessageLogger_cfi') +process.load("IOMC.EventVertexGenerators.VtxSmearedGauss_cfi") +process.load(geomFile) +process.load("Configuration.StandardSequences.MagneticField_cff") +process.load("Configuration.EventContent.EventContent_cff") +process.load('Configuration.StandardSequences.Generator_cff') +process.load('Configuration.StandardSequences.SimIdeal_cff') +process.load("SimG4CMS.Calo.CaloSimHitStudy_cfi") + +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag, globalTag, '') + +process.source = cms.Source("EmptySource") + +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(500) +) + +if 'MessageLogger' in process.__dict__: + process.MessageLogger.G4cerr=dict() + process.MessageLogger.HitStudy=dict() + +process.Timing = cms.Service("Timing") + +process.SimpleMemoryCheck = cms.Service("SimpleMemoryCheck", + oncePerEventMode = cms.untracked.bool(True), + showMallocInfo = cms.untracked.bool(True), + dump = cms.untracked.bool(True), + ignoreTotal = cms.untracked.int32(1) +) + +process.load("IOMC.RandomEngine.IOMC_cff") +process.RandomNumberGeneratorService.generator.initialSeed = 456789 +process.RandomNumberGeneratorService.g4SimHits.initialSeed = 9876 +process.RandomNumberGeneratorService.VtxSmeared.initialSeed = 123456789 +process.rndmStore = cms.EDProducer("RandomEngineStateProducer") + +process.TFileService = cms.Service("TFileService", + fileName = cms.string(histFile) +) + +# Event output +process.output = cms.OutputModule("PoolOutputModule", + process.FEVTSIMEventContent, + fileName = cms.untracked.string(outFile) +) + +process.generation_step = cms.Path(process.pgen) +process.simulation_step = cms.Path(process.psim) +process.analysis_step = cms.Path(process.CaloSimHitStudy) +process.out_step = cms.EndPath(process.output) + +process.generator.pythiaHepMCVerbosity = False +process.generator.pythiaPylistVerbosity = 0 +process.g4SimHits.Physics.type = 'SimG4Core/Physics/FTFP_BERT_EMM' +process.g4SimHits.HCalSD.TestNumberingScheme = False +process.CaloSimHitStudy.TestNumbering = False + +# Schedule definition +process.schedule = cms.Schedule(process.generation_step, + process.simulation_step, + process.analysis_step, +# process.out_step + ) + +# filter all path with the production filter sequence +for path in process.paths: + getattr(process,path)._seq = process.generator * getattr(process,path)._seq From 3645f13bef5ae258d0b9cac346740e221922122f Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sun, 6 Nov 2022 10:12:09 +0100 Subject: [PATCH 2/4] Code check --- SimG4CMS/Calo/plugins/CaloSimHitStudy.cc | 44 ++++++++++++------------ 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/SimG4CMS/Calo/plugins/CaloSimHitStudy.cc b/SimG4CMS/Calo/plugins/CaloSimHitStudy.cc index 6c7c8a7410b17..45b22519df61b 100644 --- a/SimG4CMS/Calo/plugins/CaloSimHitStudy.cc +++ b/SimG4CMS/Calo/plugins/CaloSimHitStudy.cc @@ -107,11 +107,11 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) storeRL_(ps.getUntrackedParameter("StoreRL", false)), testNumber_(ps.getUntrackedParameter("TestNumbering", true)), tokGeom_(esConsumes()), - tok_evt_(consumes(edm::InputTag(ps.getUntrackedParameter("SourceLabel", "VtxSmeared")))), - toks_calo_{edm::vector_transform(hitLab_, - [this](const std::string& name) { - return consumes(edm::InputTag{g4Label_, name}); - })} { + tok_evt_(consumes( + edm::InputTag(ps.getUntrackedParameter("SourceLabel", "VtxSmeared")))), + toks_calo_{edm::vector_transform(hitLab_, [this](const std::string& name) { + return consumes(edm::InputTag{g4Label_, name}); + })} { usesResource(TFileService::kSharedResource); for (unsigned i = 0; i != muonLab_.size(); i++) @@ -252,21 +252,21 @@ CaloSimHitStudy::CaloSimHitStudy(const edm::ParameterSet& ps) edm::LogVerbatim("HitStudy") << "CaloSimHitStudy: Completed defining histos for general tracking hits"; #endif std::string dett[nTrack_] = {"Pixel Barrel (High)", - "Pixel Endcap (High)", - "TEC (High)", - "TIB (High)", - "TID (High)", - "TOB (High)", - "Pixel Barrel (Low)", - "Pixel Endcap (Low)", - "TEC (Low)", - "TIB (Low)", - "TID (Low)", - "TOB (Low)", - "RPC", - "CSC", - "DT", - "GEM"}; + "Pixel Endcap (High)", + "TEC (High)", + "TIB (High)", + "TID (High)", + "TOB (High)", + "Pixel Barrel (Low)", + "Pixel Endcap (Low)", + "TEC (Low)", + "TIB (Low)", + "TID (Low)", + "TOB (Low)", + "RPC", + "CSC", + "DT", + "GEM"}; double nhtMax[nTrack_] = { 500., 500., 1000., 1000., 500., 1000., 5000., 2000., 10000., 5000., 2000., 5000., 500., 1000., 1000., 500.}; for (int i = 0; i < nTrack_; i++) { @@ -549,8 +549,8 @@ void CaloSimHitStudy::analyzeHits(std::vector& hits, int indx) { } if (idx > 0) { timeAll_[idx]->Fill(time); - eta_[idx]->Fill(point.eta()); - phi_[idx]->Fill(point.phi()); + eta_[idx]->Fill(point.eta()); + phi_[idx]->Fill(point.phi()); } } } From 932bc204b06eb59bd998729b0aec433017a098d1 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sun, 6 Nov 2022 12:14:24 +0100 Subject: [PATCH 3/4] Also commit the macro --- SimG4CMS/Calo/macros/MakeHitStudyPlots.C | 118 +++++++++++++---------- 1 file changed, 66 insertions(+), 52 deletions(-) diff --git a/SimG4CMS/Calo/macros/MakeHitStudyPlots.C b/SimG4CMS/Calo/macros/MakeHitStudyPlots.C index 211d55015e3ca..da69a66719f91 100644 --- a/SimG4CMS/Calo/macros/MakeHitStudyPlots.C +++ b/SimG4CMS/Calo/macros/MakeHitStudyPlots.C @@ -6,16 +6,17 @@ // 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) +// makeHitStudyPlots(file1, tag1, file2, tag2, todomin, todomax, gtitle, +// 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) +// file1 std::string Name of the first ROOT file [old/analRun3.root] +// file2 std::string Name of the second ROOT file [new/analRun3.root] +// tag1 std::string Tag for the first file [Bug] +// tag2 std::string Tag for the second file [Fix] +// gtitle std::string Overall Titile [""] +// todomin int Minimum type # of histograms to be plotted [0] +// todomax int Maximum type # of histograms to be plotted [0] // 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] @@ -42,35 +43,39 @@ #include #include -void makeHitStudyPlots(std::string file1 = "analRun3Old.root", - std::string tag1 = "Old", - std::string file2 = "analRun3New.root", - std::string tag2 = "New", - int toDo = 0, +void makeHitStudyPlots(std::string file1 = "uncorr/analRun3.root", + std::string file2 = "corr/analRun3.root", + std::string tag1 = "Bug", + std::string tag2 = "Fix", + std::string gtitle = "", + int todomin = 0, + int todomax = 0, bool ratio = true, bool save = false, std::string dirnm = "CaloSimHitStudy") { - 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}; - int todos[6] = {3, 5, 8, 9, 10, 11}; + const int plots = 20; + std::string names[plots] = {"Etot", + "Hit", + "EtotG", + "Time", + "eta", + "phi", + "EdepEM", + "EdepHad", + "EdepTk", + "Edep", + "HitHigh", + "HitLow", + "HitMu", + "HitTk", + "TimeAll", + "TimeTk", + "EneInc", + "EtaInc", + "PhiInc", + "PtInc"}; + int numb[plots] = {9, 9, 9, 9, 9, 9, 9, 9, 16, 9, 1, 1, 1, 16, 9, 16, 1, 1, 1, 1}; + int rebin[plots] = {10, 10, 10, 10, 2, 4, 10, 10, 1, 10, 10, 10, 10, 10, 10, 10, 1, 1, 1, 1}; bool debug(false); gStyle->SetCanvasBorderMode(0); @@ -100,30 +105,34 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", tag += tag2; } } - 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]; + if ((todomin < 0) || (todomin >= plots)) + todomin = 0; + if (todomax < todomin) { + todomax = todomin; + } else if (todomax >= plots) { + todomax = plots - 1; + } + std::cout << "Use " << nfile << " files from " << file1 << " and " << file2 << " and look for " << todomin << ":" + << todomax << std::endl; + for (int todo = todomin; todo <= todomax; ++todo) { 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, "c_%s%s%s", names[todo].c_str(), tag.c_str(), gtitle.c_str()); + } else { + sprintf(name, "%s%d", names[todo].c_str(), i1); + sprintf(namec, "c_%s%d%s%s", names[todo].c_str(), i1, tag.c_str(), gtitle.c_str()); + } + TCanvas* pad = new TCanvas(namec, namec, 500, 500); 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; @@ -138,6 +147,7 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", hist->GetYaxis()->SetTitleOffset(1.4); if (rebin[todo] > 1) hist->Rebin(rebin[todo]); + hist->SetTitle(gtitle.c_str()); if (first == 0) { pad = new TCanvas(namec, namec, 500, 500); pad->SetRightMargin(0.10); @@ -165,7 +175,7 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", leg->Draw("same"); pad->Update(); if (save) { - sprintf(name, "c_%s.pdf", pad->GetName()); + sprintf(name, "%s.pdf", pad->GetName()); pad->Print(name); } } @@ -176,17 +186,17 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", 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()); + sprintf(namec, "c_%sRatio%s%s", names[todo].c_str(), tag.c_str(), gtitle.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()); + sprintf(namec, "c_%s%dRatio%s%s", names[todo].c_str(), i1, tag.c_str(), gtitle.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->SetTitle(gtitle.c_str()); histr->GetXaxis()->SetTitle(hist0[0]->GetXaxis()->GetTitle()); histr->GetYaxis()->SetTitle(name); histr->GetXaxis()->SetLabelOffset(0.005); @@ -235,6 +245,10 @@ void makeHitStudyPlots(std::string file1 = "analRun3Old.root", sumDen = 0; std::cout << tag1 << " vs " << tag2 << " " << hist0[0]->GetXaxis()->GetTitle() << " Mean deviation " << sumNum << " +- " << sumDen << " maximum " << maxDev << std::endl; + if (save) { + sprintf(name, "%s.pdf", pad->GetName()); + pad->Print(name); + } } } } From e414e89f551485951d21c7d5ba013c649c890710 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sun, 6 Nov 2022 12:28:02 +0100 Subject: [PATCH 4/4] Code check --- SimG4CMS/Calo/macros/MakeHitStudyPlots.C | 45 ++++++++---------------- 1 file changed, 14 insertions(+), 31 deletions(-) diff --git a/SimG4CMS/Calo/macros/MakeHitStudyPlots.C b/SimG4CMS/Calo/macros/MakeHitStudyPlots.C index da69a66719f91..89f16d34b3527 100644 --- a/SimG4CMS/Calo/macros/MakeHitStudyPlots.C +++ b/SimG4CMS/Calo/macros/MakeHitStudyPlots.C @@ -47,33 +47,16 @@ void makeHitStudyPlots(std::string file1 = "uncorr/analRun3.root", std::string file2 = "corr/analRun3.root", std::string tag1 = "Bug", std::string tag2 = "Fix", - std::string gtitle = "", + std::string gtitle = "", int todomin = 0, int todomax = 0, bool ratio = true, bool save = false, std::string dirnm = "CaloSimHitStudy") { const int plots = 20; - std::string names[plots] = {"Etot", - "Hit", - "EtotG", - "Time", - "eta", - "phi", - "EdepEM", - "EdepHad", - "EdepTk", - "Edep", - "HitHigh", - "HitLow", - "HitMu", - "HitTk", - "TimeAll", - "TimeTk", - "EneInc", - "EtaInc", - "PhiInc", - "PtInc"}; + std::string names[plots] = {"Etot", "Hit", "EtotG", "Time", "eta", "phi", "EdepEM", + "EdepHad", "EdepTk", "Edep", "HitHigh", "HitLow", "HitMu", "HitTk", + "TimeAll", "TimeTk", "EneInc", "EtaInc", "PhiInc", "PtInc"}; int numb[plots] = {9, 9, 9, 9, 9, 9, 9, 9, 16, 9, 1, 1, 1, 16, 9, 16, 1, 1, 1, 1}; int rebin[plots] = {10, 10, 10, 10, 2, 4, 10, 10, 1, 10, 10, 10, 10, 10, 10, 10, 1, 1, 1, 1}; bool debug(false); @@ -124,11 +107,11 @@ void makeHitStudyPlots(std::string file1 = "uncorr/analRun3.root", TH1D* hist0[nfile]; char name[100], namec[100]; if (numb[todo] == 1) { - sprintf(name, "%s", names[todo].c_str()); - sprintf(namec, "c_%s%s%s", names[todo].c_str(), tag.c_str(), gtitle.c_str()); + sprintf(name, "%s", names[todo].c_str()); + sprintf(namec, "c_%s%s%s", names[todo].c_str(), tag.c_str(), gtitle.c_str()); } else { - sprintf(name, "%s%d", names[todo].c_str(), i1); - sprintf(namec, "c_%s%d%s%s", names[todo].c_str(), i1, tag.c_str(), gtitle.c_str()); + sprintf(name, "%s%d", names[todo].c_str(), i1); + sprintf(namec, "c_%s%d%s%s", names[todo].c_str(), i1, tag.c_str(), gtitle.c_str()); } TCanvas* pad = new TCanvas(namec, namec, 500, 500); for (int ifile = 0; ifile < nfile; ++ifile) { @@ -147,7 +130,7 @@ void makeHitStudyPlots(std::string file1 = "uncorr/analRun3.root", hist->GetYaxis()->SetTitleOffset(1.4); if (rebin[todo] > 1) hist->Rebin(rebin[todo]); - hist->SetTitle(gtitle.c_str()); + hist->SetTitle(gtitle.c_str()); if (first == 0) { pad = new TCanvas(namec, namec, 500, 500); pad->SetRightMargin(0.10); @@ -196,7 +179,7 @@ void makeHitStudyPlots(std::string file1 = "uncorr/analRun3.root", 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->SetTitle(gtitle.c_str()); + histr->SetTitle(gtitle.c_str()); histr->GetXaxis()->SetTitle(hist0[0]->GetXaxis()->GetTitle()); histr->GetYaxis()->SetTitle(name); histr->GetXaxis()->SetLabelOffset(0.005); @@ -245,10 +228,10 @@ void makeHitStudyPlots(std::string file1 = "uncorr/analRun3.root", sumDen = 0; std::cout << tag1 << " vs " << tag2 << " " << hist0[0]->GetXaxis()->GetTitle() << " Mean deviation " << sumNum << " +- " << sumDen << " maximum " << maxDev << std::endl; - if (save) { - sprintf(name, "%s.pdf", pad->GetName()); - pad->Print(name); - } + if (save) { + sprintf(name, "%s.pdf", pad->GetName()); + pad->Print(name); + } } } }