Merge pull request #31366 from bsunanda/Run4-hgx260

Run4-hgx260 Add a new testing code to check V14 geometry of HGCal

Validation/HGCalValidation/python/hgcalSiliconValidation_cfi.py

@@ -0,0 +1,8 @@
+import FWCore.ParameterSet.Config as cms
+
+from Validation.HGCalValidation.hgcalSiliconAnalysisEE_cfi import *
+
+hgcalSiliconAnalysisHEF = hgcalSiliconAnalysisEE.clone(
+    detectorName = cms.untracked.string("HGCalHESiliconSensitive"),
+    HitCollection = cms.untracked.string('HGCHitsHEfront'),
+    DigiCollection = cms.untracked.InputTag("hgcalDigis","HEfront"))

Validation/HGCalValidation/test/HGCalSiliconValidation.cc

@@ -0,0 +1,176 @@
+// system include files
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <map>
+#include <string>
+
+// user include files
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+#include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+#include "Geometry/Records/interface/IdealGeometryRecord.h"
+
+#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
+
+// Root objects
+#include "TROOT.h"
+#include "TSystem.h"
+#include "TFile.h"
+#include "TH1.h"
+#include "TH2.h"
+
+class HGCalSiliconValidation : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
+public:
+  explicit HGCalSiliconValidation(const edm::ParameterSet& ps);
+  ~HGCalSiliconValidation() override {}
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+protected:
+  void beginJob() override {}
+  void beginRun(edm::Run const&, edm::EventSetup const&) override;
+  void endRun(edm::Run const&, edm::EventSetup const&) override {}
+  void analyze(edm::Event const&, edm::EventSetup const&) override;
+
+private:
+  edm::Service<TFileService> fs_;
+  const std::string g4Label_, nameDetector_, hgcalHits_;
+  const edm::InputTag hgcalDigis_;
+  const int iSample_;
+  edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
+  edm::EDGetToken tok_digi_;
+
+  TH1D *hsimE1_, *hsimE2_, *hsimTm_;
+  TH1D *hsimLn_, *hdigEn_, *hdigLn_;
+  TH2D *hsimOc_, *hsi2Oc_, *hdigOc_, *hdi2Oc_;
+};
+
+HGCalSiliconValidation::HGCalSiliconValidation(const edm::ParameterSet& ps)
+    : g4Label_(ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
+      nameDetector_(ps.getUntrackedParameter<std::string>("detectorName", "HGCalEESensitive")),
+      hgcalHits_((ps.getUntrackedParameter<std::string>("HitCollection", "HGCHitsEE"))),
+      hgcalDigis_(ps.getUntrackedParameter<edm::InputTag>("DigiCollection")),
+      iSample_(ps.getUntrackedParameter<int>("Sample", 5)) {
+  usesResource(TFileService::kSharedResource);
+
+  tok_hits_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hgcalHits_));
+  tok_digi_ = consumes<HGCalDigiCollection>(hgcalDigis_);
+  edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation::Input for SimHit:"
+                                      << edm::InputTag(g4Label_, hgcalHits_) << "  Digits:" << hgcalDigis_
+                                      << "  Sample: " << iSample_;
+}
+
+void HGCalSiliconValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
+  desc.addUntracked<std::string>("detectorName", "HGCalEESensitive");
+  desc.addUntracked<std::string>("HitCollection", "HGCHitsEE");
+  desc.addUntracked<edm::InputTag>("DigiCollection", edm::InputTag("hgcalDigis", "EE"));
+  desc.addUntracked<int>("Sample", 5);
+  descriptions.add("hgcalSiliconAnalysisEE", desc);
+}
+
+void HGCalSiliconValidation::beginRun(edm::Run const&, edm::EventSetup const& es) {
+  edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation::Booking the Histograms";
+
+  //Histograms for Sim Hits
+  hsimE1_ = fs_->make<TH1D>("SimHitEn1", "Sim Hit Energy", 1000, 0.0, 1.0);
+  hsimE2_ = fs_->make<TH1D>("SimHitEn2", "Sim Hit Energy", 1000, 0.0, 1.0);
+  hsimTm_ = fs_->make<TH1D>("SimHitTime", "Sim Hit Time", 1000, 0.0, 500.0);
+  hsimLn_ = fs_->make<TH1D>("SimHitLong", "Sim Hit Long. Profile", 60, 0.0, 30.0);
+  hsimOc_ = fs_->make<TH2D>("SimHitOccup", "Sim Hit Occupnacy", 300, 0.0, 300.0, 60, 0.0, 30.0);
+  hsi2Oc_ = fs_->make<TH2D>("SimHitOccu2", "Sim Hit Occupnacy", 300, 300.0, 600.0, 300, 0.0, 300.0);
+  //Histograms for Digis
+  hdigEn_ = fs_->make<TH1D>("DigiEnergy", "Digi ADC Sample", 1000, 0.0, 1000.0);
+  hdigLn_ = fs_->make<TH1D>("DigiLong", "Digi Long. Profile", 60, 0.0, 30.0);
+  hdigOc_ = fs_->make<TH2D>("DigiOccup", "Digi Occupnacy", 300, 0.0, 300.0, 60, 0.0, 30.0);
+  hdi2Oc_ = fs_->make<TH2D>("DigiOccu2", "Digi Occupnacy", 300, 300.0, 600.0, 300, 0.0, 300.0);
+}
+
+void HGCalSiliconValidation::analyze(const edm::Event& e, const edm::EventSetup& iSetup) {
+  edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation:Run = " << e.id().run()
+                                      << " Event = " << e.id().event();
+
+  edm::ESHandle<HGCalGeometry> geom;
+  iSetup.get<IdealGeometryRecord>().get(nameDetector_, geom);
+  if (!geom.isValid()) {
+    edm::LogWarning("HGCalValidation") << "Cannot get valid HGCalGeometry Object for " << nameDetector_;
+  } else {
+    const HGCalGeometry* geom0 = geom.product();
+
+    //SimHits
+    edm::Handle<edm::PCaloHitContainer> hitsCalo;
+    e.getByToken(tok_hits_, hitsCalo);
+    edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation.: PCaloHitContainer obtained with flag "
+                                        << hitsCalo.isValid();
+    if (hitsCalo.isValid()) {
+      edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation: PCaloHit buffer " << hitsCalo->size();
+      unsigned i(0);
+      std::map<unsigned int, double> map_try;
+      for (edm::PCaloHitContainer::const_iterator it = hitsCalo->begin(); it != hitsCalo->end(); ++it) {
+        double energy = it->energy();
+        double time = it->time();
+        unsigned int id = it->id();
+        GlobalPoint pos = geom0->getPosition(DetId(id));
+        double r = pos.perp();
+        double z = std::abs(pos.z());
+        int lay = HGCSiliconDetId(id).layer();
+        hsimE1_->Fill(energy);
+        hsimTm_->Fill(time, energy);
+        hsimOc_->Fill(r, lay, energy);
+        hsi2Oc_->Fill(z, r, energy);
+        hsimLn_->Fill(lay, energy);
+        double ensum = (map_try.count(id) != 0) ? map_try[id] : 0;
+        ensum += energy;
+        map_try[id] = ensum;
+        ++i;
+        edm::LogVerbatim("HGCalValidation") << "HGCalBHHit[" << i << "] ID " << std::hex << " " << id << std::dec << " "
+                                            << HGCSiliconDetId(id) << " E " << energy << " time " << time;
+      }
+      for (std::map<unsigned int, double>::iterator itr = map_try.begin(); itr != map_try.end(); ++itr) {
+        hsimE2_->Fill((*itr).second);
+      }
+    }
+
+    //Digits
+    unsigned int kount(0);
+    edm::Handle<HGCalDigiCollection> digicoll;
+    e.getByToken(tok_digi_, digicoll);
+    edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation.: HGCalDigiCollection obtained with flag "
+                                        << digicoll.isValid();
+    if (digicoll.isValid()) {
+      edm::LogVerbatim("HGCalValidation") << "HGCalSiliconValidation: HGCalDigi buffer " << digicoll->size();
+      for (HGCalDigiCollection::const_iterator it = digicoll->begin(); it != digicoll->end(); ++it) {
+        HGCalDataFrame df(*it);
+        double energy = df[iSample_].data();
+        HGCSiliconDetId cell(df.id());
+        GlobalPoint pos = geom0->getPosition(cell);
+        double r = pos.perp();
+        double z = std::abs(pos.z());
+        int depth = cell.layer();
+        hdigEn_->Fill(energy);
+        hdigLn_->Fill(depth);
+        hdigOc_->Fill(r, depth);
+        hdi2Oc_->Fill(z, r);
+        ++kount;
+        edm::LogVerbatim("HGCalValidation") << "HGCalBHDigit[" << kount << "] ID " << cell << " E " << energy;
+      }
+    }
+  }
+}
+
+DEFINE_FWK_MODULE(HGCalSiliconValidation);

Validation/HGCalValidation/test/python/protoHGCalSimWatcher_cfg.py

@@ -0,0 +1,190 @@
+###############################################################################
+# Way to use this:
+#   cmsRun protoVHGCalSimWatcher_cfg.py geometry=D62
+#
+#   Options for geometry D49, D58, D59, D62
+#
+###############################################################################
+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',
+                 "D62",
+                  VarParsing.VarParsing.multiplicity.singleton,
+                  VarParsing.VarParsing.varType.string,
+                  "geometry of operations: D49, D58, D59, D62")
+
+### get and parse the command line arguments
+options.parseArguments()
+
+print(options)
+
+############################################################
+# Use the options
+
+if (options.geometry == "D49"):
+    from Configuration.Eras.Era_Phase2C9_cff import Phase2C9
+    process = cms.Process('PROD',Phase2C9)
+    process.load('Configuration.Geometry.GeometryExtended2026D49_cff')
+    process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff')
+elif (options.geometry == "D58"):
+    from Configuration.Eras.Era_Phase2C12_cff import Phase2C12
+    process = cms.Process('PROD',Phase2C12)
+    process.load('Configuration.Geometry.GeometryExtended2026D58_cff')
+    process.load('Configuration.Geometry.GeometryExtended2026D58Reco_cff')
+elif (options.geometry == "D59"):
+    from Configuration.Eras.Era_Phase2C11_cff import Phase2C11
+    process = cms.Process('PROD',Phase2C11)
+    process.load('Configuration.Geometry.GeometryExtended2026D59_cff')
+    process.load('Configuration.Geometry.GeometryExtended2026D59Reco_cff')
+else:
+    from Configuration.Eras.Era_Phase2C11_cff import Phase2C11
+    process = cms.Process('PROD',Phase2C11)
+    process.load('Configuration.Geometry.GeometryExtended2026D62_cff')
+    process.load('Configuration.Geometry.GeometryExtended2026D62Reco_cff')
+
+# import of standard configurations
+process.load('Configuration.StandardSequences.Services_cff')
+process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('Configuration.EventContent.EventContent_cff')
+process.load('SimGeneral.MixingModule.mixNoPU_cfi')
+process.load('Configuration.StandardSequences.MagneticField_cff')
+process.load('Configuration.StandardSequences.Generator_cff')
+process.load('IOMC.EventVertexGenerators.VtxSmearedRealistic50ns13TeVCollision_cfi')
+process.load('GeneratorInterface.Core.genFilterSummary_cff')
+process.load('Configuration.StandardSequences.SimIdeal_cff')
+process.load('Configuration.StandardSequences.EndOfProcess_cff')
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+process.load('Configuration.StandardSequences.Digi_cff')
+process.load('Configuration.StandardSequences.SimL1Emulator_cff')
+process.load('Configuration.StandardSequences.L1TrackTrigger_cff')
+process.load('Configuration.StandardSequences.DigiToRaw_cff')
+process.load('HLTrigger.Configuration.HLT_Fake2_cff')
+process.load('Configuration.StandardSequences.RawToDigi_cff')
+process.load('Configuration.StandardSequences.L1Reco_cff')
+process.load('Configuration.StandardSequences.Reconstruction_cff')
+process.load('Configuration.StandardSequences.RecoSim_cff')
+process.load('Configuration.StandardSequences.EndOfProcess_cff')
+
+process.maxEvents = cms.untracked.PSet(
+    input = cms.untracked.int32(1000)
+)
+
+process.MessageLogger.cerr.FwkReport.reportEvery = 5
+if hasattr(process,'MessageLogger'):
+    process.MessageLogger.categories.append('ValidHGCal')
+    process.MessageLogger.categories.append('HGCalGeom')
+
+# Input source
+process.source = cms.Source("EmptySource")
+
+process.options = cms.untracked.PSet(
+    wantSummary = cms.untracked.bool(True)
+)
+
+# Production Info
+process.configurationMetadata = cms.untracked.PSet(
+    version = cms.untracked.string(''),
+    annotation = cms.untracked.string(''),
+    name = cms.untracked.string('Applications')
+)
+
+# Output definition
+process.output = cms.OutputModule("PoolOutputModule",
+    splitLevel = cms.untracked.int32(0),
+    eventAutoFlushCompressedSize = cms.untracked.int32(5242880),
+    outputCommands = cms.untracked.vstring(
+        'keep *_*hbhe*_*_*',
+	'keep *_g4SimHits_*_*',
+	'keep *_*HGC*_*_*',
+        ),
+    fileName = cms.untracked.string('file:testHGCalSimWatcher.root'),
+    dataset = cms.untracked.PSet(
+        filterName = cms.untracked.string(''),
+        dataTier = cms.untracked.string('GEN-SIM-DIGI-RAW-RECO')
+    ),
+    SelectEvents = cms.untracked.PSet(
+        SelectEvents = cms.vstring('generation_step')
+    )
+)
+
+# Additional output definition
+process.g4SimHits.Watchers = cms.VPSet(cms.PSet(
+    SimG4HGCalValidation = cms.PSet(
+        Names = cms.vstring(
+                'HGCalEESensitive',  
+                'HGCalHESensitive',
+                'HGCalHESiliconSensitive',
+                'HGCalHEScintillatorSensitive',
+                ),
+        Types          = cms.vint32(0,0,0,0),
+        DetTypes       = cms.vint32(0,1,1,2),
+        LabelLayerInfo = cms.string("HGCalInfoLayer"),
+        Verbosity      = cms.untracked.int32(0),
+    ),
+    type = cms.string('SimG4HGCalValidation')
+))
+
+# Other statements
+process.genstepfilter.triggerConditions=cms.vstring("generation_step")
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:phase2_realistic', '')
+
+process.generator = cms.EDProducer("FlatRandomPtGunProducer",
+    PGunParameters = cms.PSet(
+        MaxPt = cms.double(20.0),
+        MinPt = cms.double(20.0),
+        PartID = cms.vint32(13), #--->muon
+        MaxEta = cms.double(3.0),
+        MaxPhi = cms.double(3.14159265359),
+        MinEta = cms.double(1.2),
+        MinPhi = cms.double(-3.14159265359)
+    ),
+    Verbosity = cms.untracked.int32(0),
+    psethack = cms.string('single muon pt 20'),
+    AddAntiParticle = cms.bool(False),
+    firstRun = cms.untracked.uint32(1)
+)
+
+
+#Modified to produce hgceedigis
+process.mix.digitizers = cms.PSet(process.theDigitizersValid)
+
+process.ProductionFilterSequence = cms.Sequence(process.generator)
+
+# Path and EndPath definitions
+process.generation_step = cms.Path(process.pgen)
+process.simulation_step = cms.Path(process.psim)
+process.genfiltersummary_step = cms.EndPath(process.genFilterSummary)
+process.digitisation_step = cms.Path(process.pdigi_valid)
+process.L1simulation_step = cms.Path(process.SimL1Emulator)
+process.L1TrackTrigger_step = cms.Path(process.L1TrackTrigger)
+process.digi2raw_step = cms.Path(process.DigiToRaw)
+process.raw2digi_step = cms.Path(process.RawToDigi)
+process.L1Reco_step = cms.Path(process.L1Reco)
+process.reconstruction_step = cms.Path(process.localreco)
+process.recosim_step = cms.Path(process.recosim)
+process.out_step = cms.EndPath(process.output)
+
+# Schedule definition
+process.schedule = cms.Schedule(process.generation_step,
+				process.simulation_step,
+				process.digitisation_step,
+				process.L1simulation_step,
+                                process.L1TrackTrigger_step,
+				process.digi2raw_step,
+				process.raw2digi_step,
+				process.L1Reco_step,
+				process.reconstruction_step,
+                                process.recosim_step,
+				process.out_step
+				)
+
+# filter all path with the production filter sequence
+for path in process.paths:
+        if getattr(process,path)._seq is not None: getattr(process,path)._seq = process.ProductionFilterSequence * getattr(process,path)._seq