File indexing completed on 2024-06-22 02:24:14
0001
0002 #include <cmath>
0003 #include <fstream>
0004 #include <iostream>
0005 #include <map>
0006 #include <memory>
0007 #include <string>
0008 #include <vector>
0009
0010
0011 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0012 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024
0025 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0026 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0027 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0028 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0029 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0030 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0031 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0032 #include "SimDataFormats/HcalTestBeam/interface/HcalTestBeamNumbering.h"
0033 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0034 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0035 #include "SimG4CMS/HGCalTestBeam/interface/AHCalDetId.h"
0036 #include "SimG4CMS/HGCalTestBeam/interface/AHCalGeometry.h"
0037
0038 #include "SimDataFormats/CaloHit/interface/PassiveHit.h"
0039
0040
0041 #include "TH1.h"
0042 #include "TH2.h"
0043 #include "TProfile.h"
0044 #include "TProfile2D.h"
0045 #include "TROOT.h"
0046 #include "TSystem.h"
0047 #include "TTree.h"
0048
0049
0050
0051 class HGCalTB23Analyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0052 public:
0053 explicit HGCalTB23Analyzer(edm::ParameterSet const&);
0054 ~HGCalTB23Analyzer() override = default;
0055
0056 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0057
0058 private:
0059 void beginJob() override;
0060 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0061 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0062 void analyze(edm::Event const&, edm::EventSetup const&) override;
0063 void analyzeSimHits(int type, std::vector<PCaloHit>& hits, double zFront);
0064 void analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
0065 edm::Handle<edm::SimVertexContainer> const& SimVtx);
0066 void analyzePassiveHits(edm::Handle<edm::PassiveHitContainer> const& hgcPh, int subdet);
0067 static bool sortTime(const std::pair<double, double>& i, const std::pair<double, double>& j);
0068
0069 edm::Service<TFileService> fs_;
0070 std::unique_ptr<AHCalGeometry> ahcalGeom_;
0071 const HGCalDDDConstants* hgcons_[2];
0072 const bool ifEE_, ifFH_, ifBH_, ifBeam_;
0073 const bool doSimHits_, doTree_, doTreeCell_;
0074 const bool doPassive_, doPassiveEE_, doPassiveHE_, doPassiveBH_, addP_, doBeam_;
0075 const std::string detectorEE_, detectorFH_;
0076 const std::string detectorBH_, detectorBeam_;
0077 const double zFrontEE_, zFrontFH_, zFrontBH_;
0078 const double gev2mip200_, gev2mip300_, stoc_smear_time_200_, stoc_smear_time_300_;
0079 std::vector<int> idBeams_;
0080 const edm::InputTag labelGen_;
0081 const std::string labelHitEE_, labelHitFH_, labelHitBH_, labelHitBeam_;
0082 const edm::InputTag labelPassiveEE_, labelPassiveFH_, labelPassiveBH_;
0083 const edm::InputTag labelPassiveCMSE_, labelPassiveBeam_;
0084
0085 const edm::EDGetTokenT<edm::HepMCProduct> tok_hepMC_;
0086 const edm::EDGetTokenT<edm::SimTrackContainer> tok_simTk_;
0087 const edm::EDGetTokenT<edm::SimVertexContainer> tok_simVtx_;
0088 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hitsEE_, tok_hitsFH_;
0089 const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hitsBH_, tok_hitsBeam_;
0090 const edm::EDGetTokenT<edm::PassiveHitContainer> tok_hgcPHEE_, tok_hgcPHFH_;
0091 const edm::EDGetTokenT<edm::PassiveHitContainer> tok_hgcPHBH_, tok_hgcPHCMSE_, tok_hgcPHBeam_;
0092 const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> tokDDDEE_;
0093 const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> tokDDDFH_;
0094
0095 TTree* tree_;
0096 TH1D *hSimHitE_[4], *hSimHitEn_[4], *hSimHitT_[4], *hBeam_;
0097 TProfile *hSimHitLng_[3], *hSimHitLng1_[3];
0098 TProfile* hSimHitLng2_[3];
0099 TProfile2D* hSimHitLat_[3];
0100 std::vector<TH1D*> hSimHitLayEn1EE_, hSimHitLayEn2EE_;
0101 std::vector<TH1D*> hSimHitLayEn1FH_, hSimHitLayEn2FH_;
0102 std::vector<TH1D*> hSimHitLayEn1BH_, hSimHitLayEn2BH_;
0103 std::vector<TH1D*> hSimHitLayEnBeam_;
0104 std::vector<float> simHitLayEn1EE_, simHitLayEn2EE_;
0105 std::vector<float> simHitLayEn1FH_, simHitLayEn2FH_;
0106 std::vector<float> simHitLayEn1BH_, simHitLayEn2BH_;
0107 std::vector<float> simHitLayEnBeam_;
0108 std::vector<uint32_t> simHitCellIdEE_, simHitCellIdFH_;
0109 std::vector<uint32_t> simHitCellIdBH_, simHitCellIdBeam_;
0110 std::vector<float> simHitCellEnEE_, simHitCellEnFH_;
0111 std::vector<float> simHitCellEnBH_, simHitCellEnBeam_;
0112 std::vector<int> simHitCellColBH_, simHitCellRowBH_, simHitCellLayerBH_;
0113 std::vector<float> simHitCellTimeFirstHitEE_, simHitCellTimeFirstHitFH_, simHitCellTimeFirstHitBH_;
0114 std::vector<float> simHitCellTime15MipEE_, simHitCellTime15MipFH_, simHitCellTime15MipBH_;
0115 std::vector<float> simHitCellTimeLastHitEE_, simHitCellTimeLastHitFH_, simHitCellTimeLastHitBH_;
0116
0117 std::vector<float> hgcPassiveEEEnergy_, hgcPassiveFHEnergy_, hgcPassiveBHEnergy_, hgcPassiveCMSEEnergy_;
0118 std::vector<float> hgcPassiveBeamEnergy_;
0119 std::vector<std::string> hgcPassiveEEName_, hgcPassiveFHName_, hgcPassiveBHName_, hgcPassiveCMSEName_;
0120 std::vector<std::string> hgcPassiveBeamName_;
0121 std::vector<int> hgcPassiveEEID_, hgcPassiveFHID_, hgcPassiveBHID_, hgcPassiveCMSEID_, hgcPassiveBeamID_;
0122
0123 double xBeam_, yBeam_, zBeam_, pBeam_;
0124 double thetaBeam_, phiBeam_;
0125 int nBeamMC_;
0126 std::vector<int> pdgIdBeamMC_;
0127 std::vector<float> xBeamMC_, yBeamMC_, zBeamMC_;
0128 std::vector<float> pxBeamMC_, pyBeamMC_, pzBeamMC_, pBeamMC_;
0129 };
0130
0131 HGCalTB23Analyzer::HGCalTB23Analyzer(const edm::ParameterSet& iConfig)
0132 : ifEE_(iConfig.getParameter<bool>("useEE")),
0133 ifFH_(iConfig.getParameter<bool>("useFH")),
0134 ifBH_(iConfig.getParameter<bool>("useBH")),
0135 ifBeam_(iConfig.getParameter<bool>("useBeam")),
0136 doSimHits_(iConfig.getParameter<bool>("doSimHits")),
0137 doTree_(iConfig.getParameter<bool>("doTree")),
0138 doTreeCell_(iConfig.getParameter<bool>("doTreeCell")),
0139 doPassive_(iConfig.getParameter<bool>("doPassive")),
0140 doPassiveEE_(iConfig.getParameter<bool>("doPassiveEE")),
0141 doPassiveHE_(iConfig.getParameter<bool>("doPassiveHE")),
0142 doPassiveBH_(iConfig.getParameter<bool>("doPassiveBH")),
0143 addP_(iConfig.getParameter<bool>("addP")),
0144 doBeam_(iConfig.getParameter<bool>("doBeam")),
0145 detectorEE_(iConfig.getParameter<std::string>("detectorEE")),
0146 detectorFH_(iConfig.getParameter<std::string>("detectorFH")),
0147 detectorBH_(iConfig.getParameter<std::string>("detectorBH")),
0148 detectorBeam_(iConfig.getParameter<std::string>("detectorBeam")),
0149 zFrontEE_(iConfig.getParameter<double>("zFrontEE")),
0150 zFrontFH_(iConfig.getParameter<double>("zFrontFH")),
0151 zFrontBH_(iConfig.getParameter<double>("zFrontBH")),
0152 gev2mip200_(iConfig.getUntrackedParameter<double>("gev2mip200", 57.0e-6)),
0153 gev2mip300_(iConfig.getUntrackedParameter<double>("gev2mip300", 85.5e-6)),
0154 stoc_smear_time_200_(iConfig.getUntrackedParameter<double>("stoc_smear_time_200", 10.24)),
0155 stoc_smear_time_300_(iConfig.getUntrackedParameter<double>("stoc_smear_time_300", 15.5)),
0156 labelGen_(iConfig.getParameter<edm::InputTag>("generatorSrc")),
0157 labelHitEE_(iConfig.getParameter<std::string>("caloHitSrcEE")),
0158 labelHitFH_(iConfig.getParameter<std::string>("caloHitSrcFH")),
0159 labelHitBH_(iConfig.getParameter<std::string>("caloHitSrcBH")),
0160 labelHitBeam_(iConfig.getParameter<std::string>("caloHitSrcBeam")),
0161 labelPassiveEE_(iConfig.getParameter<edm::InputTag>("passiveEE")),
0162 labelPassiveFH_(iConfig.getParameter<edm::InputTag>("passiveFH")),
0163 labelPassiveBH_(iConfig.getParameter<edm::InputTag>("passiveBH")),
0164 labelPassiveCMSE_(iConfig.getParameter<edm::InputTag>("passiveCMSE")),
0165 labelPassiveBeam_(iConfig.getParameter<edm::InputTag>("passiveBeam")),
0166 tok_hepMC_(consumes<edm::HepMCProduct>(labelGen_)),
0167 tok_simTk_(consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"))),
0168 tok_simVtx_(consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"))),
0169 tok_hitsEE_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitEE_))),
0170 tok_hitsFH_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitFH_))),
0171 tok_hitsBH_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitBH_))),
0172 tok_hitsBeam_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", labelHitBeam_))),
0173 tok_hgcPHEE_(consumes<edm::PassiveHitContainer>(labelPassiveEE_)),
0174 tok_hgcPHFH_(consumes<edm::PassiveHitContainer>(labelPassiveFH_)),
0175 tok_hgcPHBH_(consumes<edm::PassiveHitContainer>(labelPassiveBH_)),
0176 tok_hgcPHCMSE_(consumes<edm::PassiveHitContainer>(labelPassiveCMSE_)),
0177 tok_hgcPHBeam_(consumes<edm::PassiveHitContainer>(labelPassiveBeam_)),
0178 tokDDDEE_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0179 edm::ESInputTag("", detectorEE_))),
0180 tokDDDFH_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0181 edm::ESInputTag("", detectorFH_))) {
0182 usesResource("TFileService");
0183 ahcalGeom_ = std::make_unique<AHCalGeometry>(iConfig);
0184
0185
0186
0187
0188
0189
0190
0191 idBeams_ = (iConfig.getParameter<std::vector<int>>("idBeams"));
0192 #ifdef EDM_ML_DEBUG
0193 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: SimHits = " << doSimHits_ << " useDets " << ifEE_ << ":" << ifFH_
0194 << ":" << ifBH_ << ":" << ifBeam_ << " zFront " << zFrontEE_ << ":" << zFrontFH_ << ":"
0195 << zFrontBH_ << " IdBeam " << idBeams_.size() << ":";
0196 for (unsigned int k = 0; k < idBeams_.size(); ++k)
0197 edm::LogVerbatim("HGCSim") << " [" << k << "] " << idBeams_[k];
0198 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: DoPassive " << doPassive_ << ":" << doPassiveEE_ << ":"
0199 << doPassiveHE_ << ":" << doPassiveBH_;
0200 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: MIP conversion factors " << gev2mip200_ << ":" << gev2mip300_
0201 << " Time smearing " << stoc_smear_time_200_ << ":" << stoc_smear_time_300_ << " AddP "
0202 << addP_;
0203 #endif
0204 if (idBeams_.empty())
0205 idBeams_.push_back(1001);
0206
0207 #ifdef EDM_ML_DEBUG
0208 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: GeneratorSource = " << labelGen_;
0209 if (ifEE_)
0210 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: Detector " << detectorEE_ << " with tags " << labelHitEE_;
0211 if (ifFH_)
0212 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: Detector " << detectorFH_ << " with tags " << labelHitFH_;
0213 if (ifBH_)
0214 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: Detector " << detectorBH_ << " with tags " << labelHitBH_;
0215 if (ifBeam_)
0216 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: Detector " << detectorBeam_ << " with tags " << labelHitBeam_;
0217 #endif
0218 }
0219
0220 void HGCalTB23Analyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0221 edm::ParameterSetDescription desc;
0222 desc.add<std::string>("detectorEE", "HGCalEESensitive");
0223 desc.add<bool>("useEE", true);
0224 desc.add<double>("zFrontEE", 0.0);
0225 desc.add<std::string>("caloHitSrcEE", "HGCHitsEE");
0226 desc.add<std::string>("detectorFH", "HGCalHESiliconSensitive");
0227 desc.add<bool>("useFH", false);
0228 desc.add<double>("zFrontFH", 0.0);
0229 desc.add<std::string>("caloHitSrcFH", "HGCHitsHEfront");
0230 desc.add<std::string>("detectorBH", "AHCal");
0231 desc.add<bool>("useBH", false);
0232 desc.add<double>("zFrontBH", 0.0);
0233 desc.add<std::string>("caloHitSrcBH", "HcalHits");
0234 desc.add<std::string>("detectorBeam", "HcalTB06BeamDetector");
0235 desc.add<bool>("useBeam", false);
0236 desc.add<std::string>("caloHitSrcBeam", "HcalTB06BeamHits");
0237 std::vector<int> ids = {
0238 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1011, 1012, 1013, 1014, 2001, 2002, 2003, 2004, 2005};
0239 desc.add<std::vector<int>>("idBeams", ids);
0240 desc.add<edm::InputTag>("generatorSrc", edm::InputTag("generatorSmeared"));
0241 desc.add<edm::InputTag>("passiveEE", edm::InputTag("g4SimHits", "HGCalEEPassiveHits"));
0242 desc.add<edm::InputTag>("passiveFH", edm::InputTag("g4SimHits", "HGCalHEPassiveHits"));
0243 desc.add<edm::InputTag>("passiveBH", edm::InputTag("g4SimHits", "HGCalAHPassiveHits"));
0244 desc.add<edm::InputTag>("passiveCMSE", edm::InputTag("g4SimHits", "CMSEPassiveHits"));
0245 desc.add<edm::InputTag>("passiveBeam", edm::InputTag("g4SimHits", "HGCalBeamPassiveHits"));
0246
0247 desc.add<bool>("doSimHits", true);
0248 desc.add<bool>("doTree", true);
0249 desc.add<bool>("doTreeCell", true);
0250 desc.add<bool>("doPassive", false);
0251 desc.add<bool>("doPassiveEE", false);
0252 desc.add<bool>("doPassiveHE", false);
0253 desc.add<bool>("doPassiveBH", false);
0254 desc.add<bool>("addP", false);
0255 desc.add<bool>("doBeam", false);
0256 desc.addUntracked<double>("gev2mip200", 57.0e-6);
0257 desc.addUntracked<double>("gev2mip300", 85.5e-6);
0258 desc.addUntracked<double>("stoc_smear_time_200", 10.24);
0259 desc.addUntracked<double>("stoc_smear_time_300", 15.5);
0260 desc.addUntracked<int>("maxDepth", 12);
0261 desc.addUntracked<double>("deltaX", 30.0);
0262 desc.addUntracked<double>("deltaY", 30.0);
0263 desc.addUntracked<double>("deltaZ", 81.0);
0264 desc.addUntracked<double>("zFirst", 17.6);
0265 descriptions.add("HGCalTB23Analyzer", desc);
0266 }
0267
0268 void HGCalTB23Analyzer::beginJob() {
0269 char name[40], title[100];
0270 hBeam_ = fs_->make<TH1D>("BeamP", "Beam Momentum", 1000, 0, 1000.0);
0271 for (int i = 0; i < 3; ++i) {
0272 bool book(ifEE_);
0273 std::string det(detectorEE_);
0274 if (i == 1) {
0275 book = ifFH_;
0276 det = detectorFH_;
0277 } else if (i == 2) {
0278 book = ifBH_;
0279 det = detectorBH_;
0280 }
0281
0282 if (doSimHits_ && book) {
0283 sprintf(name, "SimHitEn%s", det.c_str());
0284 sprintf(title, "Sim Hit Energy for %s", det.c_str());
0285 hSimHitE_[i] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
0286 sprintf(name, "SimHitEnX%s", det.c_str());
0287 sprintf(title, "Sim Hit Energy for %s", det.c_str());
0288 hSimHitEn_[i] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
0289 sprintf(name, "SimHitTm%s", det.c_str());
0290 sprintf(title, "Sim Hit Timing for %s", det.c_str());
0291 hSimHitT_[i] = fs_->make<TH1D>(name, title, 5000, 0., 500.0);
0292 sprintf(name, "SimHitLat%s", det.c_str());
0293 sprintf(title, "Lateral Shower profile (Sim Hit) for %s", det.c_str());
0294 hSimHitLat_[i] = fs_->make<TProfile2D>(name, title, 100, -100., 100., 100, -100., 100.);
0295 sprintf(name, "SimHitLng%s", det.c_str());
0296 sprintf(title, "Longitudinal Shower profile (Sim Hit) for %s", det.c_str());
0297 hSimHitLng_[i] = fs_->make<TProfile>(name, title, 50, 0., 100.);
0298 sprintf(name, "SimHitLng1%s", det.c_str());
0299 sprintf(title, "Longitudinal Shower profile (Layer) for %s", det.c_str());
0300 hSimHitLng1_[i] = fs_->make<TProfile>(name, title, 200, 0., 100.);
0301 sprintf(name, "SimHitLng2%s", det.c_str());
0302 sprintf(title, "Longitudinal Shower profile (Layer) for %s", det.c_str());
0303 hSimHitLng2_[i] = fs_->make<TProfile>(name, title, 200, 0., 100.);
0304 }
0305 }
0306 if (ifBeam_ && doSimHits_) {
0307 sprintf(name, "SimHitEn%s", detectorBeam_.c_str());
0308 sprintf(title, "Sim Hit Energy for %s", detectorBeam_.c_str());
0309 hSimHitE_[3] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
0310 sprintf(name, "SimHitEnX%s", detectorBeam_.c_str());
0311 sprintf(title, "Sim Hit Energy for %s", detectorBeam_.c_str());
0312 hSimHitEn_[3] = fs_->make<TH1D>(name, title, 100000, 0., 0.2);
0313 sprintf(name, "SimHitTm%s", detectorBeam_.c_str());
0314 sprintf(title, "Sim Hit Timing for %s", detectorBeam_.c_str());
0315 hSimHitT_[3] = fs_->make<TH1D>(name, title, 5000, 0., 500.0);
0316 }
0317 if (doSimHits_ && doTree_) {
0318 tree_ = fs_->make<TTree>("HGCTB", "SimHitEnergy");
0319 tree_->Branch("simHitLayEn1EE", &simHitLayEn1EE_);
0320 tree_->Branch("simHitLayEn2EE", &simHitLayEn2EE_);
0321 tree_->Branch("simHitLayEn1FH", &simHitLayEn1FH_);
0322 tree_->Branch("simHitLayEn2FH", &simHitLayEn2FH_);
0323 tree_->Branch("simHitLayEn1BH", &simHitLayEn1BH_);
0324 tree_->Branch("simHitLayEn2BH", &simHitLayEn2BH_);
0325 tree_->Branch("xBeam", &xBeam_, "xBeam/D");
0326 tree_->Branch("yBeam", &yBeam_, "yBeam/D");
0327 tree_->Branch("zBeam", &zBeam_, "zBeam/D");
0328 tree_->Branch("pBeam", &pBeam_, "pBeam/D");
0329 tree_->Branch("thetaBeam", &thetaBeam_, "thetaBeam/D");
0330 tree_->Branch("phiBeam", &phiBeam_, "phiBeam/D");
0331 if (doBeam_) {
0332 tree_->Branch("nBeamMC", &nBeamMC_, "nBeamMC/I");
0333 tree_->Branch("pdgIdBeamMC", &pdgIdBeamMC_);
0334 tree_->Branch("xBeamMC", &xBeamMC_);
0335 tree_->Branch("yBeamMC", &yBeamMC_);
0336 tree_->Branch("zBeamMC", &zBeamMC_);
0337 tree_->Branch("pxBeamMC", &pxBeamMC_);
0338 tree_->Branch("pyBeamMC", &pyBeamMC_);
0339 tree_->Branch("pzBeamMC", &pzBeamMC_);
0340 tree_->Branch("pBeamMC", &pBeamMC_);
0341 }
0342 if (doTreeCell_) {
0343 tree_->Branch("simHitCellIdEE", &simHitCellIdEE_);
0344 tree_->Branch("simHitCellEnEE", &simHitCellEnEE_);
0345 tree_->Branch("simHitCellIdFH", &simHitCellIdFH_);
0346 tree_->Branch("simHitCellEnFH", &simHitCellEnFH_);
0347 tree_->Branch("simHitCellIdBH", &simHitCellIdBH_);
0348 tree_->Branch("simHitCellEnBH", &simHitCellEnBH_);
0349 tree_->Branch("simHitCellIdBeam", &simHitCellIdBeam_);
0350 tree_->Branch("simHitCellEnBeam", &simHitCellEnBeam_);
0351
0352 tree_->Branch("simHitCellColBH", &simHitCellColBH_);
0353 tree_->Branch("simHitCellRowBH", &simHitCellRowBH_);
0354 tree_->Branch("simHitCellLayerBH", &simHitCellLayerBH_);
0355 tree_->Branch("simHitCellTimeFirstHitEE", &simHitCellTimeFirstHitEE_);
0356 tree_->Branch("simHitCellTimeFirstHitFH", &simHitCellTimeFirstHitFH_);
0357
0358 tree_->Branch("simHitCellTime15MipEE", &simHitCellTime15MipEE_);
0359 tree_->Branch("simHitCellTime15MipFH", &simHitCellTime15MipFH_);
0360
0361 tree_->Branch("simHitCellTimeLastHitEE", &simHitCellTimeLastHitEE_);
0362 tree_->Branch("simHitCellTimeLastHitFH", &simHitCellTimeLastHitFH_);
0363
0364 }
0365 }
0366
0367 if (doPassive_ && doTree_) {
0368 if (doPassiveEE_) {
0369 tree_->Branch("hgcPassiveEEEnergy", &hgcPassiveEEEnergy_);
0370 tree_->Branch("hgcPassiveEEName", &hgcPassiveEEName_);
0371 tree_->Branch("hgcPassiveEEID", &hgcPassiveEEID_);
0372 }
0373 if (doPassiveHE_) {
0374 tree_->Branch("hgcPassiveFHEnergy", &hgcPassiveFHEnergy_);
0375 tree_->Branch("hgcPassiveFHName", &hgcPassiveFHName_);
0376 tree_->Branch("hgcPassiveFHID", &hgcPassiveFHID_);
0377 }
0378 if (doPassiveBH_) {
0379 tree_->Branch("hgcPassiveBHEnergy", &hgcPassiveBHEnergy_);
0380 tree_->Branch("hgcPassiveBHName", &hgcPassiveBHName_);
0381 tree_->Branch("hgcPassiveBHID", &hgcPassiveBHID_);
0382 }
0383 tree_->Branch("hgcPassiveCMSEEnergy", &hgcPassiveCMSEEnergy_);
0384 tree_->Branch("hgcPassiveCMSEName", &hgcPassiveCMSEName_);
0385 tree_->Branch("hgcPassiveCMSEID", &hgcPassiveCMSEID_);
0386 tree_->Branch("hgcPassiveBeamEnergy", &hgcPassiveBeamEnergy_);
0387 tree_->Branch("hgcPassiveBeamName", &hgcPassiveBeamName_);
0388 tree_->Branch("hgcPassiveBeamID", &hgcPassiveBeamID_);
0389 }
0390 }
0391
0392 void HGCalTB23Analyzer::beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0393 char name[40], title[100];
0394 if (ifEE_) {
0395 hgcons_[0] = &iSetup.getData(tokDDDEE_);
0396 for (unsigned int l = 0; l < hgcons_[0]->layers(false); ++l) {
0397 sprintf(name, "SimHitEnA%d%s", l, detectorEE_.c_str());
0398 sprintf(title, "Sim Hit Energy in SIM layer %d for %s", l + 1, detectorEE_.c_str());
0399 hSimHitLayEn1EE_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0400 if (l % 3 == 0) {
0401 sprintf(name, "SimHitEnB%d%s", (l / 3 + 1), detectorEE_.c_str());
0402 sprintf(title, "Sim Hit Energy in layer %d for %s", (l / 3 + 1), detectorEE_.c_str());
0403 hSimHitLayEn2EE_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0404 }
0405 }
0406 #ifdef EDM_ML_DEBUG
0407 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer::" << detectorEE_ << " defined with " << hgcons_[0]->layers(false)
0408 << " layers";
0409 #endif
0410 } else {
0411 hgcons_[0] = nullptr;
0412 }
0413
0414 if (ifFH_) {
0415 hgcons_[1] = &iSetup.getData(tokDDDFH_);
0416 for (unsigned int l = 0; l < hgcons_[1]->layers(false); ++l) {
0417 sprintf(name, "SimHitEnA%d%s", l, detectorFH_.c_str());
0418 sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorFH_.c_str());
0419 hSimHitLayEn1FH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0420 if (l % 3 == 0) {
0421 sprintf(name, "SimHitEnB%d%s", (l / 3 + 1), detectorFH_.c_str());
0422 sprintf(title, "Sim Hit Energy in layer %d for %s", (l / 3 + 1), detectorFH_.c_str());
0423 hSimHitLayEn2FH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0424 }
0425 }
0426 #ifdef EDM_ML_DEBUG
0427 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer::" << detectorFH_ << " defined with " << hgcons_[1]->layers(false)
0428 << " layers";
0429 #endif
0430 } else {
0431 hgcons_[1] = nullptr;
0432 }
0433
0434 if (ifBH_) {
0435 for (int l = 0; l < ahcalGeom_->maxDepth(); ++l) {
0436 sprintf(name, "SimHitEnA%d%s", l, detectorBH_.c_str());
0437 sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorBH_.c_str());
0438 hSimHitLayEn1BH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0439 sprintf(name, "SimHitEnB%d%s", l, detectorBH_.c_str());
0440 sprintf(title, "Sim Hit Energy in layer %d for %s", l + 1, detectorBH_.c_str());
0441 hSimHitLayEn2BH_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0442 }
0443 }
0444
0445 if (ifBeam_) {
0446 for (unsigned int l = 0; l < idBeams_.size(); ++l) {
0447 sprintf(name, "SimHitEna%d%s", l, detectorBeam_.c_str());
0448 sprintf(title, "Sim Hit Energy in type %d for %s", idBeams_[l], detectorBeam_.c_str());
0449 hSimHitLayEnBeam_.push_back(fs_->make<TH1D>(name, title, 100000, 0., 0.2));
0450 }
0451 }
0452 }
0453
0454 void HGCalTB23Analyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0455
0456 const edm::Handle<edm::HepMCProduct>& evtMC = iEvent.getHandle(tok_hepMC_);
0457 if (!evtMC.isValid()) {
0458 edm::LogWarning("HGCal") << "no HepMCProduct found";
0459 } else {
0460 const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
0461 unsigned int k(0);
0462 HepMC::FourVector pxyz(0, 0, 0, 0);
0463 for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0464 ++p, ++k) {
0465 edm::LogVerbatim("HGCSim") << "Particle [" << k << "] with p " << (*p)->momentum().rho() << " theta "
0466 << (*p)->momentum().theta() << " phi " << (*p)->momentum().phi() << " pxyz ("
0467 << (*p)->momentum().px() << ", " << (*p)->momentum().py() << ", "
0468 << (*p)->momentum().pz() << ")";
0469 if (addP_) {
0470 pxyz.setPx(pxyz.px() + (*p)->momentum().px());
0471 pxyz.setPy(pxyz.py() + (*p)->momentum().py());
0472 pxyz.setPz(pxyz.pz() + (*p)->momentum().pz());
0473 pxyz.setE(pxyz.e() + (*p)->momentum().e());
0474 } else if (!addP_ && (k == 0)) {
0475 pxyz = (*p)->momentum();
0476 }
0477 }
0478 hBeam_->Fill(pxyz.rho());
0479 edm::LogVerbatim("HGCSim") << "Particle with p " << pxyz.rho() << " theta " << pxyz.theta() << " phi "
0480 << pxyz.phi();
0481 }
0482
0483
0484 if (doSimHits_) {
0485 const edm::Handle<edm::SimTrackContainer>& SimTk = iEvent.getHandle(tok_simTk_);
0486 const edm::Handle<edm::SimVertexContainer>& SimVtx = iEvent.getHandle(tok_simVtx_);
0487 analyzeSimTracks(SimTk, SimVtx);
0488
0489 simHitLayEn1EE_.clear();
0490 simHitLayEn2EE_.clear();
0491 simHitLayEn1FH_.clear();
0492 simHitLayEn2FH_.clear();
0493 simHitLayEn1BH_.clear();
0494 simHitLayEn2BH_.clear();
0495 simHitLayEnBeam_.clear();
0496 simHitCellIdEE_.clear();
0497 simHitCellEnEE_.clear();
0498 simHitCellIdFH_.clear();
0499 simHitCellEnFH_.clear();
0500 simHitCellIdBH_.clear();
0501 simHitCellEnBH_.clear();
0502 simHitCellIdBeam_.clear();
0503 simHitCellEnBeam_.clear();
0504 simHitCellColBH_.clear();
0505 simHitCellRowBH_.clear();
0506 simHitCellLayerBH_.clear();
0507 simHitCellTimeFirstHitEE_.clear();
0508 simHitCellTime15MipEE_.clear();
0509 simHitCellTimeLastHitEE_.clear();
0510 simHitCellTimeFirstHitFH_.clear();
0511 simHitCellTime15MipFH_.clear();
0512 simHitCellTimeLastHitFH_.clear();
0513 std::vector<PCaloHit> caloHits;
0514 if (ifEE_) {
0515 simHitLayEn1EE_ = std::vector<float>(hgcons_[0]->layers(false), 0);
0516 simHitLayEn2EE_ = std::vector<float>(hgcons_[0]->layers(true), 0);
0517 const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsEE_);
0518 if (theCaloHitContainers.isValid()) {
0519 #ifdef EDM_ML_DEBUG
0520 edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorEE_ << " has " << theCaloHitContainers->size()
0521 << " hits";
0522 #endif
0523 caloHits.clear();
0524 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0525 analyzeSimHits(0, caloHits, zFrontEE_);
0526 } else {
0527 #ifdef EDM_ML_DEBUG
0528 edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorEE_ << " !!!";
0529 #endif
0530 }
0531 }
0532 if (ifFH_) {
0533 simHitLayEn1FH_ = std::vector<float>(hgcons_[1]->layers(false), 0);
0534 simHitLayEn2FH_ = std::vector<float>(hgcons_[1]->layers(true), 0);
0535 const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsFH_);
0536 if (theCaloHitContainers.isValid()) {
0537 #ifdef EDM_ML_DEBUG
0538 edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorFH_ << " has " << theCaloHitContainers->size()
0539 << " hits";
0540 #endif
0541 caloHits.clear();
0542 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0543 analyzeSimHits(1, caloHits, zFrontFH_);
0544 } else {
0545 #ifdef EDM_ML_DEBUG
0546 edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorFH_ << " !!!";
0547 #endif
0548 }
0549 }
0550 if (ifBH_) {
0551 simHitLayEn1BH_ = std::vector<float>(ahcalGeom_->maxDepth(), 0);
0552 simHitLayEn2BH_ = std::vector<float>(ahcalGeom_->maxDepth(), 0);
0553 const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsBH_);
0554 if (theCaloHitContainers.isValid()) {
0555 #ifdef EDM_ML_DEBUG
0556 edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorBH_ << " has " << theCaloHitContainers->size()
0557 << " hits";
0558 #endif
0559 caloHits.clear();
0560 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0561 analyzeSimHits(2, caloHits, zFrontBH_);
0562 } else {
0563 #ifdef EDM_ML_DEBUG
0564 edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorBH_ << " !!!";
0565 #endif
0566 }
0567 }
0568 if (ifBeam_) {
0569 simHitLayEnBeam_ = std::vector<float>(idBeams_.size(), 0);
0570 const edm::Handle<edm::PCaloHitContainer>& theCaloHitContainers = iEvent.getHandle(tok_hitsBeam_);
0571 if (theCaloHitContainers.isValid()) {
0572 #ifdef EDM_ML_DEBUG
0573 edm::LogVerbatim("HGCSim") << "PcalohitContainer for " << detectorBeam_ << " has "
0574 << theCaloHitContainers->size() << " hits";
0575 #endif
0576 caloHits.clear();
0577 caloHits.insert(caloHits.end(), theCaloHitContainers->begin(), theCaloHitContainers->end());
0578 analyzeSimHits(3, caloHits, 0.0);
0579 } else {
0580 #ifdef EDM_ML_DEBUG
0581 edm::LogVerbatim("HGCSim") << "PCaloHitContainer does not exist for " << detectorBeam_ << " !!!";
0582 #endif
0583 }
0584 }
0585 }
0586
0587
0588 if (doPassive_) {
0589
0590 hgcPassiveEEEnergy_.clear();
0591 hgcPassiveEEName_.clear();
0592 hgcPassiveEEID_.clear();
0593 if (doPassiveEE_) {
0594 const edm::Handle<edm::PassiveHitContainer>& hgcPHEE = iEvent.getHandle(tok_hgcPHEE_);
0595 analyzePassiveHits(hgcPHEE, 1);
0596 }
0597
0598 hgcPassiveFHEnergy_.clear();
0599 hgcPassiveFHName_.clear();
0600 hgcPassiveFHID_.clear();
0601 if (doPassiveHE_) {
0602 const edm::Handle<edm::PassiveHitContainer>& hgcPHFH = iEvent.getHandle(tok_hgcPHFH_);
0603 analyzePassiveHits(hgcPHFH, 2);
0604 }
0605
0606 hgcPassiveBHEnergy_.clear();
0607 hgcPassiveBHName_.clear();
0608 hgcPassiveBHID_.clear();
0609 if (doPassiveBH_) {
0610 const edm::Handle<edm::PassiveHitContainer>& hgcPHBH = iEvent.getHandle(tok_hgcPHBH_);
0611 analyzePassiveHits(hgcPHBH, 3);
0612 }
0613
0614 hgcPassiveCMSEEnergy_.clear();
0615 hgcPassiveCMSEName_.clear();
0616 hgcPassiveCMSEID_.clear();
0617 const edm::Handle<edm::PassiveHitContainer>& hgcPHCMSE = iEvent.getHandle(tok_hgcPHCMSE_);
0618 analyzePassiveHits(hgcPHCMSE, 4);
0619
0620 hgcPassiveBeamName_.clear();
0621 hgcPassiveBeamEnergy_.clear();
0622 hgcPassiveBeamID_.clear();
0623 const edm::Handle<edm::PassiveHitContainer>& hgcPHBeam = iEvent.getHandle(tok_hgcPHBeam_);
0624 analyzePassiveHits(hgcPHBeam, 5);
0625 }
0626
0627 if ((doSimHits_ || doPassive_) && (doTree_))
0628 tree_->Fill();
0629 }
0630
0631 void HGCalTB23Analyzer::analyzeSimHits(int type, std::vector<PCaloHit>& hits, double zFront) {
0632 std::map<uint32_t, double> map_hits, map_hitn;
0633 std::map<uint32_t, double> map_hittime_firsthit, map_hittime_lasthit, map_hittime_15Mip;
0634 std::map<int, double> map_hitDepth, map_hitWafer;
0635 std::map<int, std::pair<uint32_t, double>> map_hitLayer, map_hitCell;
0636 double entot(0);
0637 std::map<uint32_t, double> nhits;
0638 std::map<uint32_t, int> ID, Depth;
0639 std::map<uint32_t, double> GeV2Mip;
0640 std::map<uint32_t, double> StochTermTime;
0641 std::vector<int> nSimLayers;
0642
0643 std::map<uint32_t, std::vector<std::pair<double, double>>> map_hitTimeEn;
0644
0645 bool debug = false;
0646 for (unsigned int i = 0; i < hits.size(); i++) {
0647 double energy = hits[i].energy();
0648 double time = hits[i].time();
0649 uint32_t id = hits[i].id();
0650 entot += energy;
0651 int subdet, zside, layer, sector, subsector(0), cell, cell2(0), depth(0), idx(0);
0652 if (type == 2) {
0653 subdet = HcalDetId(id).subdet();
0654 if (subdet != HcalOther)
0655 continue;
0656 AHCalDetId hid(id);
0657 layer = depth = hid.depth();
0658 zside = hid.zside();
0659 sector = hid.irow();
0660 cell = hid.icol();
0661 idx = ((hid.irowAbs() * 100) + (hid.icolAbs()));
0662 if (debug)
0663 edm::LogVerbatim("HGCSim") << "depth, sector, cell " << depth << ":" << sector << ":" << cell;
0664 } else if (type == 3) {
0665 HcalTestBeamNumbering::unpackIndex(id, subdet, layer, sector, cell);
0666 depth = layer;
0667 zside = 1;
0668 idx = subdet * 1000 + layer;
0669 layer = idx;
0670 } else {
0671 #ifdef EDM_ML_DEBUG
0672 edm::LogVerbatim("HGCSim") << "ID " << std::hex << id << std::dec << " " << HGCSiliconDetId(id);
0673 #endif
0674 HGCSiliconDetId(id).unpack(subdet, zside, layer, sector, subsector, cell, cell2);
0675 #ifdef EDM_ML_DEBUG
0676 edm::LogVerbatim("HGCSim") << "Unpack ID " << subdet << " " << zside << " " << layer << " " << sector << " "
0677 << subsector << " " << cell << " " << cell2;
0678 #endif
0679 depth = hgcons_[type]->simToReco(cell, layer, sector, true).second;
0680 static constexpr uint32_t mask = 0xFFFFF;
0681 static constexpr uint32_t maskcell = 0x3FF;
0682 idx = (id & mask);
0683 #ifdef EDM_ML_DEBUG
0684 std::pair<float, float> xy = hgcons_[type]->locateCell(HGCSiliconDetId(id), false);
0685 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer::detId " << std::hex << id << std::dec << " Layer:Wafer:Cell "
0686 << layer << ":" << sector << ":" << subsector << ":" << cell << ":" << cell2
0687 << " Position " << xy.first << ":" << xy.second << ":"
0688 << hgcons_[type]->waferZ(layer, false);
0689 #endif
0690 cell = (idx & maskcell);
0691 sector = (idx >> 10) & maskcell;
0692 }
0693 #ifdef EDM_ML_DEBUG
0694 edm::LogVerbatim("HGCSim") << "SimHit:Hit[" << i << "] Id " << subdet << ":" << zside << ":" << layer << ":"
0695 << sector << ":" << cell << " Energy " << energy << " Time " << time;
0696 #endif
0697 if (map_hits.count(id) != 0) {
0698 map_hits[id] += energy;
0699 } else {
0700 map_hits[id] = energy;
0701 }
0702 if (map_hitLayer.count(layer) != 0) {
0703 double ee = energy + map_hitLayer[layer].second;
0704 map_hitLayer[layer] = std::make_pair(id, ee);
0705 } else {
0706 map_hitLayer[layer] = std::make_pair(id, energy);
0707 }
0708 if (map_hitWafer.count(sector) != 0)
0709 map_hitWafer[sector] += energy;
0710 else
0711 map_hitWafer[sector] = energy;
0712 if (depth >= 0) {
0713 if (map_hitCell.count(idx) != 0) {
0714 double ee = energy + map_hitCell[idx].second;
0715 map_hitCell[idx] = std::make_pair(id, ee);
0716 } else {
0717 map_hitCell[idx] = std::make_pair(id, energy);
0718 }
0719 if (debug) {
0720 if (type == 1)
0721 edm::LogVerbatim("HGCSim") << "EE, depth is and map_hitDepth[depth] " << depth << " " << map_hitDepth[depth];
0722 if (type == 2)
0723 edm::LogVerbatim("HGCSim") << "BH, depth is and map_hitDepth[depth] " << depth << " " << map_hitDepth[depth];
0724 }
0725 if (map_hitDepth.count(depth) != 0) {
0726 map_hitDepth[depth] += energy;
0727 } else {
0728 map_hitDepth[depth] = energy;
0729 }
0730 uint32_t idn =
0731 (type >= 2) ? id : HGCalTestNumbering::packHexagonIndex(subdet, zside, depth, sector, subsector, cell);
0732 map_hitTimeEn[idn].push_back(std::make_pair(time, energy));
0733 GeV2Mip[idn] = gev2mip300_;
0734 StochTermTime[idn] = stoc_smear_time_300_;
0735 ID[idn] = id;
0736 Depth[idn] = depth;
0737 if (map_hitn.count(idn) != 0) {
0738 map_hitn[idn] += energy;
0739 ++nhits[idn];
0740 } else {
0741 map_hitn[idn] = energy;
0742 nhits[idn] = 1;
0743 }
0744 }
0745 hSimHitT_[type]->Fill(time, energy);
0746 }
0747
0748 if (type < 2) {
0749 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer:: " << map_hitWafer.size() << " wafers are hit in type " << type;
0750 for (auto itr = map_hitWafer.begin(); itr != map_hitWafer.end(); ++itr)
0751 edm::LogVerbatim("HGCSim") << "Wafer: " << itr->first << " Deposited Energy " << itr->second;
0752
0753 for (const auto& itr : map_hitTimeEn) {
0754 uint32_t id = itr.first;
0755 int waferU(-99), waferV(-99);
0756
0757 waferU = HGCSiliconDetId(ID[id]).waferU();
0758 waferV = HGCSiliconDetId(ID[id]).waferV();
0759 double layer = HGCSiliconDetId(id).layer();
0760 double thickness = hgcons_[type]->cellThickness(layer, waferU, waferV);
0761 if (debug)
0762 edm::LogVerbatim("HGCSim") << "wafer is : depth (reco) " << waferU << ":" << waferV << " " << Depth[id]
0763 << "\ntype : layer : wafer thickness " << type << " " << layer << " " << thickness
0764 << "\nID(sim) and id(reco) " << std::hex << ID[id] << " " << id << std::dec;
0765 if (thickness == 300) {
0766 GeV2Mip[id] = gev2mip300_;
0767 StochTermTime[id] = stoc_smear_time_300_;
0768 } else if (thickness == 200) {
0769 GeV2Mip[id] = gev2mip200_;
0770 StochTermTime[id] = stoc_smear_time_200_;
0771 }
0772
0773 std::sort(map_hitTimeEn[id].begin(), map_hitTimeEn[id].end(), sortTime);
0774
0775 double threshold = 15.;
0776 double totE = 0.;
0777 double totEbeforeThreshold = 0.;
0778 double timebeforeThreshold = 0.;
0779 double timeAtThresohld = 0.;
0780 for (unsigned int ihit = 0; ihit < map_hitTimeEn[id].size(); ihit++) {
0781 double energy = (map_hitTimeEn[id].at(ihit)).second / GeV2Mip[id];
0782 totE += energy;
0783 double time = (map_hitTimeEn[id].at(ihit)).first;
0784 if (debug)
0785 edm::LogVerbatim("HGCSim") << "Tot E till now : time of that E : GeV2Mip[id] is " << totE << " " << time
0786 << " " << GeV2Mip[id];
0787 if (totE < threshold) {
0788 totEbeforeThreshold = totE;
0789 timebeforeThreshold = time;
0790 } else {
0791 timeAtThresohld =
0792 (threshold - totEbeforeThreshold) * (time - timebeforeThreshold) / (totE - totEbeforeThreshold) +
0793 timebeforeThreshold;
0794 map_hittime_15Mip[id] = timeAtThresohld;
0795 if (debug)
0796 edm::LogVerbatim("HGCSim") << "ihit : energyBefore : timeBefore : energyTot : timeTot : timeAt15MIP "
0797 << ihit << " " << totEbeforeThreshold << " " << timebeforeThreshold << " "
0798 << totE << " " << time << " " << map_hittime_15Mip[id];
0799 break;
0800 }
0801 }
0802 if (!map_hitTimeEn[id].empty()) {
0803 map_hittime_firsthit[id] = (map_hitTimeEn[id].at(0)).first;
0804 map_hittime_lasthit[id] = (map_hitTimeEn[id].at(map_hitTimeEn[id].size() - 1)).first;
0805 if (map_hittime_15Mip[id] < map_hittime_firsthit[id])
0806 map_hittime_15Mip[id] = map_hittime_firsthit[id];
0807
0808
0809
0810
0811
0812
0813 }
0814
0815 if (totE < threshold)
0816 map_hittime_15Mip[id] = -99;
0817 if (debug)
0818 edm::LogVerbatim("HGCSim") << "id : first hit time : last hit time " << id << " " << map_hittime_firsthit[id]
0819 << " " << map_hittime_lasthit[id] << "\nFinally for this cell, time is "
0820 << map_hittime_15Mip[id];
0821 }
0822 }
0823
0824 hSimHitEn_[type]->Fill(entot);
0825 for (const auto& itr : map_hits) {
0826 hSimHitE_[type]->Fill(itr.second);
0827 }
0828
0829 if (debug)
0830 edm::LogVerbatim("HGCSim") << "Now looping over map_hitLayer";
0831 for (const auto& itr : map_hitLayer) {
0832 int layer = (type == 2) ? itr.first : (itr.first - 1);
0833 double energy = (itr.second).second;
0834 double zp(0);
0835 if (type < 2)
0836 zp = hgcons_[type]->waferZ(layer + 1, false);
0837 else if (type == 2)
0838 zp = ahcalGeom_->getZ(AHCalDetId((itr.second).first));
0839 #ifdef EDM_ML_DEBUG
0840 edm::LogVerbatim("HGCSim") << "SimHit:Layer " << layer + 1 << " Z " << zp << ":" << zp - zFront << " E " << energy;
0841 #endif
0842 if (type < 3) {
0843 hSimHitLng_[type]->Fill(zp - zFront, energy);
0844 hSimHitLng2_[type]->Fill(layer + 1, energy);
0845 }
0846 if (type == 0) {
0847 if (layer < static_cast<int>(hSimHitLayEn1EE_.size())) {
0848 simHitLayEn1EE_[layer] = energy;
0849 hSimHitLayEn1EE_[layer]->Fill(energy);
0850 }
0851 } else if (type == 1) {
0852 if (layer < static_cast<int>(hSimHitLayEn1FH_.size())) {
0853 simHitLayEn1FH_[layer] = energy;
0854 hSimHitLayEn1FH_[layer]->Fill(energy);
0855 }
0856 } else if (type == 2) {
0857 if (debug)
0858 edm::LogVerbatim("HGCSim") << "layer < hSimHitLayEn1BH_.size()";
0859 if (layer < static_cast<int>(hSimHitLayEn1BH_.size())) {
0860 simHitLayEn1BH_[layer] = energy;
0861 hSimHitLayEn1BH_[layer]->Fill(energy);
0862 }
0863 } else {
0864 for (unsigned int k = 0; k < idBeams_.size(); ++k) {
0865 if (layer + 1 == idBeams_[k]) {
0866 simHitLayEnBeam_[k] = energy;
0867 hSimHitLayEnBeam_[k]->Fill(energy);
0868 break;
0869 }
0870 }
0871 }
0872 }
0873 for (const auto& itr : map_hitDepth) {
0874 int layer = (type == 2) ? itr.first : (itr.first - 1);
0875 double energy = itr.second;
0876 #ifdef EDM_ML_DEBUG
0877 edm::LogVerbatim("HGCSim") << "SimHit:Layer " << layer + 1 << " " << energy;
0878 #endif
0879 hSimHitLng1_[type]->Fill(layer + 1, energy);
0880 if (type == 0) {
0881 if (layer < static_cast<int>(hSimHitLayEn2EE_.size())) {
0882 simHitLayEn2EE_[layer] = energy;
0883 hSimHitLayEn2EE_[layer]->Fill(energy);
0884 }
0885 } else if (type == 1) {
0886 if (layer < static_cast<int>(hSimHitLayEn2FH_.size())) {
0887 simHitLayEn2FH_[layer] = energy;
0888 hSimHitLayEn2FH_[layer]->Fill(energy);
0889 }
0890 } else if (type == 2) {
0891 if (debug)
0892 edm::LogVerbatim("HGCSim") << "Inside map_hitDepth, layer no. " << layer;
0893 if (layer < static_cast<int>(hSimHitLayEn2BH_.size())) {
0894 simHitLayEn2BH_[layer] = energy;
0895 hSimHitLayEn2BH_[layer]->Fill(energy);
0896 }
0897 }
0898 }
0899
0900 if (debug)
0901 edm::LogVerbatim("HGCSim") << "Now looping over map_hitCell";
0902 if (type < 3) {
0903 for (const auto& itr : map_hitCell) {
0904 uint32_t id = ((itr.second).first);
0905 double energy = ((itr.second).second);
0906 std::pair<float, float> xy(0, 0);
0907 double xx(0);
0908 if (type == 2) {
0909 xy = ahcalGeom_->getXY(AHCalDetId(id));
0910 xx = xy.first;
0911 } else {
0912 int subdet, zside, layer, sector, subsector, cell;
0913 HGCalTestNumbering::unpackHexagonIndex(id, subdet, zside, layer, sector, subsector, cell);
0914 xy = hgcons_[type]->locateCell(cell, layer, sector, false);
0915 double zp = hgcons_[type]->waferZ(layer, false);
0916 xx = (zp < 0) ? -xy.first : xy.first;
0917 }
0918 hSimHitLat_[type]->Fill(xx, xy.second, energy);
0919 }
0920 }
0921
0922 for (const auto& itr : map_hitn) {
0923 uint32_t id = itr.first;
0924 double energy = itr.second;
0925 if (type == 0) {
0926 double time_firsthit = map_hittime_firsthit[id];
0927 double time15Mip = map_hittime_15Mip[id];
0928 double time_lasthit = map_hittime_lasthit[id];
0929 simHitCellIdEE_.push_back(id);
0930 simHitCellEnEE_.push_back(energy);
0931 simHitCellTimeFirstHitEE_.push_back(time_firsthit);
0932 simHitCellTime15MipEE_.push_back(time15Mip);
0933 simHitCellTimeLastHitEE_.push_back(time_lasthit);
0934 if (debug && (energy / GeV2Mip[id] < 15) && (map_hittime_15Mip[id] > 0))
0935 edm::LogVerbatim("HGCSim") << "FOUND!!!!rechit energy : Finally for this cell, time is " << energy / GeV2Mip[id]
0936 << " " << map_hittime_15Mip[id];
0937 } else if (type == 1) {
0938 double time_firsthit = map_hittime_firsthit[id];
0939 double time15Mip = map_hittime_15Mip[id];
0940 double time_lasthit = map_hittime_lasthit[id];
0941 simHitCellIdFH_.push_back(id);
0942 simHitCellEnFH_.push_back(energy);
0943 simHitCellTimeFirstHitFH_.push_back(time_firsthit);
0944 simHitCellTime15MipFH_.push_back(time15Mip);
0945 simHitCellTimeLastHitFH_.push_back(time_lasthit);
0946 } else if (type == 2) {
0947 simHitCellIdBH_.push_back(id);
0948 simHitCellEnBH_.push_back(energy);
0949 AHCalDetId hid(id);
0950 int row = hid.irow();
0951 int col = hid.icol();
0952 int layer = hid.depth();
0953 simHitCellColBH_.push_back(col);
0954 simHitCellRowBH_.push_back(row);
0955 simHitCellLayerBH_.push_back(layer);
0956 #ifdef EDM_ML_DEBUG
0957 edm::LogVerbatim("HGCSim") << "ID: " << std::hex << id << std::dec << " Layer: " << layer << " col: " << col
0958 << " row: " << row;
0959 #endif
0960 } else if (type == 3) {
0961 simHitCellIdBeam_.push_back(id);
0962 simHitCellEnBeam_.push_back(energy);
0963 }
0964 }
0965 }
0966
0967 void HGCalTB23Analyzer::analyzeSimTracks(edm::Handle<edm::SimTrackContainer> const& SimTk,
0968 edm::Handle<edm::SimVertexContainer> const& SimVtx) {
0969 xBeam_ = yBeam_ = zBeam_ = pBeam_ = -9999;
0970 nBeamMC_ = thetaBeam_ = phiBeam_ = -9999;
0971 int nParBeam = 0;
0972 int vertIndex(-1);
0973 if (doBeam_) {
0974 pdgIdBeamMC_.clear();
0975 xBeamMC_.clear();
0976 yBeamMC_.clear();
0977 zBeamMC_.clear();
0978 pxBeamMC_.clear();
0979 pyBeamMC_.clear();
0980 pzBeamMC_.clear();
0981 pBeamMC_.clear();
0982 }
0983 std::vector<float> verX, verY, verZ;
0984 verX.clear();
0985 verY.clear();
0986 verZ.clear();
0987 for (const auto& simVtxItr : *SimVtx) {
0988 verX.push_back(simVtxItr.position().X());
0989 verY.push_back(simVtxItr.position().Y());
0990 verZ.push_back(simVtxItr.position().Z());
0991 }
0992 #ifdef EDM_ML_DEBUG
0993 edm::LogVerbatim("HGCSim") << "Size of track " << SimTk->size();
0994 #endif
0995 HepMC::FourVector pxyz(0, 0, 0, 0);
0996 for (const auto& simTrkItr : *SimTk) {
0997 if (addP_ && !(simTrkItr.noGenpart())) {
0998 pxyz.setPx(pxyz.px() + simTrkItr.momentum().px());
0999 pxyz.setPy(pxyz.py() + simTrkItr.momentum().py());
1000 pxyz.setPz(pxyz.pz() + simTrkItr.momentum().pz());
1001 pxyz.setE(pxyz.e() + simTrkItr.momentum().e());
1002 #ifdef EDM_ML_DEBUG
1003 edm::LogVerbatim("HGCSim") << "Track " << simTrkItr.trackId() << " Vertex " << simTrkItr.vertIndex() << " Type "
1004 << simTrkItr.type() << " Charge " << simTrkItr.charge() << " px "
1005 << simTrkItr.momentum().px() << " py " << simTrkItr.momentum().py() << " pz "
1006 << simTrkItr.momentum().pz() << " P " << simTrkItr.momentum().P() << " GenIndex "
1007 << simTrkItr.genpartIndex();
1008 edm::LogVerbatim("HGCSim") << "Vertex " << simTrkItr.vertIndex()
1009 << " position-> X: " << verX[simTrkItr.vertIndex()]
1010 << " Y: " << verY[simTrkItr.vertIndex()] << " Z: " << verZ[simTrkItr.vertIndex()];
1011 #endif
1012 }
1013 if (doBeam_ && !(simTrkItr.noGenpart())) {
1014 nParBeam++;
1015 pdgIdBeamMC_.push_back(simTrkItr.type());
1016 xBeamMC_.push_back(verX[simTrkItr.vertIndex()]);
1017 yBeamMC_.push_back(verY[simTrkItr.vertIndex()]);
1018 zBeamMC_.push_back(verZ[simTrkItr.vertIndex()]);
1019 pxBeamMC_.push_back(simTrkItr.momentum().px());
1020 pyBeamMC_.push_back(simTrkItr.momentum().py());
1021 pzBeamMC_.push_back(simTrkItr.momentum().pz());
1022 pBeamMC_.push_back(simTrkItr.momentum().P());
1023 } else if (!addP_ && (vertIndex == -1)) {
1024 pxyz = simTrkItr.momentum();
1025 }
1026 if (vertIndex == -1)
1027 vertIndex = simTrkItr.vertIndex();
1028 }
1029 nBeamMC_ = nParBeam;
1030 pBeam_ = pxyz.rho();
1031 thetaBeam_ = pxyz.theta();
1032 phiBeam_ = pxyz.phi();
1033 if (phiBeam_ < 0)
1034 phiBeam_ += (2 * M_PI);
1035 if (vertIndex != -1 && vertIndex < static_cast<int>(SimVtx->size())) {
1036 edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
1037 for (int iv = 0; iv < vertIndex; iv++)
1038 simVtxItr++;
1039 edm::LogVerbatim("HGCSim") << "Vertex " << vertIndex << " position " << simVtxItr->position();
1040 xBeam_ = verX[0];
1041 yBeam_ = verY[0];
1042 zBeam_ = verZ[0];
1043 }
1044 }
1045
1046 void HGCalTB23Analyzer::analyzePassiveHits(edm::Handle<edm::PassiveHitContainer> const& hgcPH, int subdet) {
1047 for (const auto& v : *hgcPH) {
1048 double energy = v.energy();
1049 std::string name = v.vname();
1050 unsigned int id = v.id();
1051 #ifdef EDM_ML_DEBUG
1052 double time = v.time();
1053 edm::LogVerbatim("HGCSim") << "HGCalTB23Analyzer::analyzePassiveHits:Energy:"
1054 << "Time:Name:Id : " << energy << ":" << time << ":" << name << ":" << id;
1055 #endif
1056
1057 if (subdet == 1) {
1058 hgcPassiveEEEnergy_.push_back(energy);
1059 hgcPassiveEEName_.push_back(name);
1060 hgcPassiveEEID_.push_back(id);
1061 } else if (subdet == 2) {
1062 hgcPassiveFHEnergy_.push_back(energy);
1063 hgcPassiveFHName_.push_back(name);
1064 hgcPassiveFHID_.push_back(id);
1065 } else if (subdet == 3) {
1066 hgcPassiveBHEnergy_.push_back(energy);
1067 hgcPassiveBHName_.push_back(name);
1068 hgcPassiveBHID_.push_back(id);
1069 } else if (subdet == 4) {
1070 hgcPassiveCMSEEnergy_.push_back(energy);
1071 hgcPassiveCMSEName_.push_back(name);
1072 hgcPassiveCMSEID_.push_back(id);
1073 } else if (subdet == 5) {
1074 hgcPassiveBeamEnergy_.push_back(energy);
1075 hgcPassiveBeamName_.push_back(name);
1076 hgcPassiveBeamID_.push_back(id);
1077 }
1078 }
1079 }
1080
1081 bool HGCalTB23Analyzer::sortTime(const std::pair<double, double>& i, const std::pair<double, double>& j) {
1082 return i.first < j.first;
1083 }
1084
1085
1086 DEFINE_FWK_MODULE(HGCalTB23Analyzer);
