HcalSimAlgos/test/HcalDigitizerTest.cc

0001 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0002 #include "DataFormats/HcalDigi/interface/HBHEDataFrame.h"
0003 #include "DataFormats/HcalDigi/interface/HODataFrame.h"
0004 #include "DataFormats/HcalDigi/interface/HFDataFrame.h"
0005 #include "DataFormats/HcalDigi/interface/ZDCDataFrame.h"
0006 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0007 #include "SimCalorimetry/CaloSimAlgos/interface/CaloHitResponse.h"
0008 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSiPMHitResponse.h"
0009 #include "SimCalorimetry/CaloSimAlgos/interface/CaloTDigitizer.h"
0010 #include "SimCalorimetry/CaloSimAlgos/interface/CaloShapeIntegrator.h"
0011 #include "SimCalorimetry/CaloSimAlgos/interface/CaloShapes.h"
0012 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameterMap.h"
0013 #include "SimCalorimetry/HcalSimAlgos/interface/HcalShape.h"
0014 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSiPMShape.h"
0015 #include "SimCalorimetry/HcalSimAlgos/interface/HFShape.h"
0016 #include "SimCalorimetry/HcalSimAlgos/interface/ZDCShape.h"
0017 #include "SimCalorimetry/HcalSimAlgos/interface/HcalElectronicsSim.h"
0018 #include "CalibCalorimetry/HcalAlgos/interface/HcalDbHardcode.h"
0019 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0020 #include "SimCalorimetry/HcalSimAlgos/interface/HcalAmplifier.h"
0021 #include "SimCalorimetry/HcalSimAlgos/interface/HcalCoderFactory.h"
0022 #include "SimCalorimetry/HcalSimAlgos/interface/HcalHitFilter.h"
0023 #include "SimCalorimetry/HcalSimAlgos/interface/ZDCHitFilter.h"
0024 #include "CondFormats/DataRecord/interface/HcalTimeSlewRecord.h"
0025 #include "CondFormats/HcalObjects/interface/HcalPedestals.h"
0026 #include "CondFormats/HcalObjects/interface/HcalPedestalWidths.h"
0027 #include "CondFormats/HcalObjects/interface/HcalGains.h"
0028 #include "CondFormats/HcalObjects/interface/HcalGainWidths.h"
0029 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0030 #include "SimCalorimetry/HcalSimAlgos/interface/HcalDigitizerTraits.h"
0031 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0032 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0033 #include "Geometry/HcalTowerAlgo/interface/HcalHardcodeGeometryLoader.h"
0034 #include "FWCore/Framework/interface/Frameworkfwd.h"
0035 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0039 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0040 #include "CLHEP/Random/JamesRandom.h"
0041 #include <vector>
0042 #include <iostream>
0043 #include <iterator>
0044
0045 class HcalTimeSlew;
0046
0047 class HcalDigitizerTest : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0048 public:
0049   explicit HcalDigitizerTest(const edm::ParameterSet& iConfig);
0050   ~HcalDigitizerTest() override;
0051
0052 private:
0053   void beginJob() override;
0054   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0055   void analyze(edm::Event const&, edm::EventSetup const&) override;
0056   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0057
0058   HcalDbHardcode dbHardcode;
0059   std::vector<PCaloHit> hits;
0060   std::vector<DetId> hcalDetIds, hoDetIds, hfDetIds, hzdcDetIds, allDetIds;
0061   std::vector<HcalDetId> outerHcalDetIds;
0062
0063   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0064   const edm::ESGetToken<HcalTimeSlew, HcalTimeSlewRecord> tok_slew_;
0065
0066   const HcalTimeSlew* hcalTimeSlew_delay_;
0067 };
0068
0069 HcalDigitizerTest::HcalDigitizerTest(const edm::ParameterSet& iConfig)
0070     : tok_htopo_(esConsumes<HcalTopology, HcalRecNumberingRecord>()),
0071       tok_slew_(esConsumes<HcalTimeSlew, HcalTimeSlewRecord>(edm::ESInputTag("", "HBHE"))) {
0072   //DB helper preparation
0073   dbHardcode.setHB(HcalHardcodeParameters(iConfig.getParameter<edm::ParameterSet>("hb")));
0074   dbHardcode.setHE(HcalHardcodeParameters(iConfig.getParameter<edm::ParameterSet>("he")));
0075   dbHardcode.setHF(HcalHardcodeParameters(iConfig.getParameter<edm::ParameterSet>("hf")));
0076   dbHardcode.setHO(HcalHardcodeParameters(iConfig.getParameter<edm::ParameterSet>("ho")));
0077   dbHardcode.setHBUpgrade(HcalHardcodeParameters(iConfig.getParameter<edm::ParameterSet>("hbUpgrade")));
0078   dbHardcode.setHEUpgrade(HcalHardcodeParameters(iConfig.getParameter<edm::ParameterSet>("heUpgrade")));
0079   dbHardcode.setHFUpgrade(HcalHardcodeParameters(iConfig.getParameter<edm::ParameterSet>("hfUpgrade")));
0080   dbHardcode.useHBUpgrade(iConfig.getParameter<bool>("useHBUpgrade"));
0081   dbHardcode.useHEUpgrade(iConfig.getParameter<bool>("useHEUpgrade"));
0082   dbHardcode.useHFUpgrade(iConfig.getParameter<bool>("useHFUpgrade"));
0083   dbHardcode.testHFQIE10(iConfig.getParameter<bool>("testHFQIE10"));
0084
0085   hcalTimeSlew_delay_ = nullptr;
0086 }
0087
0088 HcalDigitizerTest::~HcalDigitizerTest() {}
0089
0090 void HcalDigitizerTest::beginJob() {
0091   // make a silly little hit in each subdetector, which should
0092   // correspond to a 100 GeV particle
0093
0094   for (int phi = 1; phi < 50; ++phi) {
0095     HcalDetId detId(HcalBarrel, 1, phi, 1);
0096     PCaloHit barrelHit(detId.rawId(), 0.085 * phi, 0.);
0097     hcalDetIds.push_back(detId);
0098     hits.push_back(barrelHit);
0099   }
0100
0101   HcalDetId endcapDetId(HcalEndcap, 17, 1, 1);
0102   PCaloHit endcapHit(endcapDetId.rawId(), 0.9, 0.);
0103   hcalDetIds.push_back(endcapDetId);
0104   hits.push_back(endcapHit);
0105
0106   HcalDetId outerDetId(HcalOuter, 1, 1, 4);
0107   PCaloHit outerHit(outerDetId.rawId(), 0.45, 0.);
0108   hoDetIds.push_back(outerDetId);
0109   outerHcalDetIds.push_back(outerDetId);
0110   hits.push_back(outerHit);
0111
0112   HcalDetId forwardDetId1(HcalForward, 30, 1, 1);
0113   PCaloHit forwardHit1(forwardDetId1.rawId(), 35.2, 0.);
0114   hfDetIds.push_back(forwardDetId1);
0115   hits.push_back(forwardHit1);
0116
0117   HcalDetId forwardDetId2(HcalForward, 30, 1, 2);
0118   PCaloHit forwardHit2(forwardDetId2.rawId(), 47.8, 0.);
0119   hfDetIds.push_back(forwardDetId2);
0120   hits.push_back(forwardHit2);
0121
0122   HcalZDCDetId zdcDetId(HcalZDCDetId::Section(1), true, 1);
0123   PCaloHit zdcHit(zdcDetId.rawId(), 50.0, 0.123);
0124   hzdcDetIds.push_back(zdcDetId);
0125   hits.push_back(zdcHit);
0126
0127   std::cout << zdcDetId << std::endl;
0128   std::cout << zdcHit << std::endl;
0129
0130   allDetIds.insert(allDetIds.end(), hcalDetIds.begin(), hcalDetIds.end());
0131   allDetIds.insert(allDetIds.end(), hoDetIds.begin(), hoDetIds.end());
0132   allDetIds.insert(allDetIds.end(), hfDetIds.begin(), hfDetIds.end());
0133   allDetIds.insert(allDetIds.end(), hzdcDetIds.begin(), hzdcDetIds.end());
0134 }
0135
0136 void HcalDigitizerTest::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0137   auto topology = iSetup.getData(tok_htopo_);
0138   hcalTimeSlew_delay_ = &iSetup.getData(tok_slew_);
0139
0140   std::string hitsName = "HcalHits";
0141   std::vector<std::string> caloDets;
0142   CrossingFrame<PCaloHit> crossingFrame(-5, 5, 25, hitsName, 0);
0143   crossingFrame.addSignals(&hits, iEvent.id());
0144
0145   // make 1 GeV pileup hit
0146   HcalDetId barrelDetId(HcalBarrel, 1, 1, 1);
0147   PCaloHit barrelPileup(barrelDetId.rawId(), 0.00855, 0.);
0148   // 10 GeV pileup hit
0149   HcalDetId forwardDetId1(HcalForward, 30, 1, 1);
0150   PCaloHit forwardPileup(forwardDetId1.rawId(), 3.52, 0.);
0151   HcalZDCDetId zdcDetId(HcalZDCDetId::Section(1), true, 1);
0152   PCaloHit zdcPileup(zdcDetId.rawId(), 3.52, 0.);
0153
0154   std::vector<PCaloHit> pileups;
0155   pileups.push_back(barrelPileup);
0156   pileups.push_back(forwardPileup);
0157   pileups.push_back(zdcPileup);
0158   ///TODO fix once the new crossingframe is released
0159   //crossingFrame.addPileupCaloHits(-3, hitsName, &pileups);
0160   // -or -
0161   // crossingFrame.addPileupCaloHits(-3, hitsName, &pileups, 0);
0162   HcalSimParameterMap parameterMap;
0163   HcalSimParameterMap siPMParameterMap = parameterMap;
0164   siPMParameterMap.setHOZecotekDetIds(outerHcalDetIds);
0165   HcalShape hcalShape;
0166   HcalSiPMShape sipmShape;
0167   HFShape hfShape;
0168   ZDCShape zdcShape;
0169
0170   CaloShapeIntegrator hcalShapeIntegrator(new HcalShape());
0171   CaloShapeIntegrator sipmShapeIntegrator(new HcalSiPMShape());
0172   CaloShapeIntegrator hfShapeIntegrator(new HFShape());
0173   CaloShapeIntegrator zdcShapeIntegrator(new ZDCShape());
0174   CaloShapes sipmShapes(&sipmShapeIntegrator);
0175   //for(float t = -25; t < 200; t += 5) {
0176   //  std::cout <<  t << " " << hcalShape(t) << "  " << sipmShape(t) << "  " << hcalShapeIntegrator(t) << "  "<< sipmShapeIntegrator(t) << std::endl;
0177   //}
0178
0179   CaloHitResponse hbheResponse(&parameterMap, &hcalShapeIntegrator);
0180   CaloHitResponse hoResponse(&parameterMap, &hcalShapeIntegrator);
0181   CaloHitResponse hfResponse(&parameterMap, &hfShapeIntegrator);
0182   CaloHitResponse zdcResponse(&parameterMap, &zdcShapeIntegrator);
0183   HcalSiPMHitResponse hoSiPMResponse(&siPMParameterMap, &sipmShapes);
0184
0185   CLHEP::HepJamesRandom randomEngine;
0186
0187   HBHEHitFilter hbheHitFilter;
0188   HOHitFilter hoHitFilter;
0189   HFHitFilter hfHitFilter;
0190   ZDCHitFilter zdcHitFilter;
0191
0192   hbheResponse.setHitFilter(&hbheHitFilter);
0193   hoSiPMResponse.setHitFilter(&hoHitFilter);
0194   hoResponse.setHitFilter(&hoHitFilter);
0195   hfResponse.setHitFilter(&hfHitFilter);
0196   zdcResponse.setHitFilter(&zdcHitFilter);
0197
0198   HcalPedestals pedestals(&topology);
0199   HcalPedestalWidths pedestalWidths(&topology);
0200   HcalGains gains(&topology);
0201   HcalGainWidths gainWidths(&topology);
0202   // make a calibration service by hand
0203   for (auto detItr = allDetIds.begin(); detItr != allDetIds.end(); ++detItr) {
0204     pedestals.addValues(dbHardcode.makePedestal(*detItr, false, false, NULL, 0.0));
0205     pedestalWidths.addValues(dbHardcode.makePedestalWidth(*detItr, false, NULL, 0.0));
0206     gains.addValues(dbHardcode.makeGain(*detItr));
0207     gainWidths.addValues(dbHardcode.makeGainWidth(*detItr));
0208   }
0209
0210   //pedestals.sort();
0211   //pedestalWidths.sort();
0212   //gains.sort();
0213   //gainWidths.sort();
0214
0215   //std::cout << "TEST Pedestal " << pedestals.getValues(barrelDetId,  1) << std::endl;
0216   //std::cout << "ZDC pedestal " << pedestals.getValue(zdcDetId,  1) << std::endl;
0217   //std::cout << "ZDC pedestal width " << pedestalWidths.getWidth(zdcDetId,  1) << std::endl;
0218   //std::cout << "ZDC gain " << gains.getValue(zdcDetId,  1) << std::endl;
0219   //std::cout << "ZDC gain width " << gainWidths.getValue(zdcDetId,  1) << std::endl;
0220
0221   HcalDbService calibratorHandle;
0222   calibratorHandle.setData(&pedestals);
0223   calibratorHandle.setData(&pedestalWidths);
0224   calibratorHandle.setData(&gains);
0225   calibratorHandle.setData(&gainWidths);
0226
0227   bool addNoise = false;
0228   bool PM1 = false;
0229   bool PM2 = false;
0230   HcalAmplifier amplifier(&parameterMap, addNoise, PM1, PM2);
0231   HcalCoderFactory coderFactory(HcalCoderFactory::NOMINAL);
0232   HcalElectronicsSim electronicsSim(&parameterMap, &amplifier, &coderFactory, PM1);
0233   amplifier.setDbService(&calibratorHandle);
0234   amplifier.setTimeSlew(hcalTimeSlew_delay_);
0235   parameterMap.setDbService(&calibratorHandle);
0236   siPMParameterMap.setDbService(&calibratorHandle);
0237
0238   CaloTDigitizer<HBHEDigitizerTraits> hbheDigitizer(&hbheResponse, &electronicsSim, addNoise);
0239   CaloTDigitizer<HODigitizerTraits> hoDigitizer(&hoResponse, &electronicsSim, addNoise);
0240   CaloTDigitizer<HFDigitizerTraits> hfDigitizer(&hfResponse, &electronicsSim, addNoise);
0241   CaloTDigitizer<ZDCDigitizerTraits> zdcDigitizer(&zdcResponse, &electronicsSim, addNoise);
0242
0243   hbheDigitizer.setDetIds(hcalDetIds);
0244   hfDigitizer.setDetIds(hfDetIds);
0245   hoDigitizer.setDetIds(hoDetIds);
0246   zdcDigitizer.setDetIds(hzdcDetIds);
0247
0248   std::unique_ptr<HBHEDigiCollection> hbheResult(new HBHEDigiCollection);
0249   std::unique_ptr<HODigiCollection> hoResult(new HODigiCollection);
0250   std::unique_ptr<HFDigiCollection> hfResult(new HFDigiCollection);
0251   std::unique_ptr<ZDCDigiCollection> zdcResult(new ZDCDigiCollection);
0252
0253   MixCollection<PCaloHit> hitCollection(&crossingFrame);
0254
0255   std::cout << "HBHE " << std::endl;
0256   hbheResponse.run(hitCollection, &randomEngine);
0257   std::cout << "SIPM " << std::endl;
0258   hoSiPMResponse.run(hitCollection, &randomEngine);
0259   //hbheDigitizer.run(hitCollection, *hbheResult);
0260   //hoDigitizer.run(hitCollection, *hoResult);
0261   //hfDigitizer.run(hitCollection, *hfResult);
0262   //zdcDigitizer.run(hitCollection, *zdcResult);
0263
0264   // print out all the digis
0265   std::cout << "HBHE Frames" << std::endl;
0266   std::copy(hbheResult->begin(), hbheResult->end(), std::ostream_iterator<HBHEDataFrame>(std::cout, "\n"));
0267
0268   std::cout << "HF Frames" << std::endl;
0269   std::copy(hfResult->begin(), hfResult->end(), std::ostream_iterator<HFDataFrame>(std::cout, "\n"));
0270
0271   std::cout << "ZDC Frames" << std::endl;
0272   std::copy(zdcResult->begin(), zdcResult->end(), std::ostream_iterator<ZDCDataFrame>(std::cout, "\n"));
0273 }
0274
0275 //define this as a plug-in
0276 DEFINE_FWK_MODULE(HcalDigitizerTest);