File indexing completed on 2023-03-17 11:24:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <iostream>
0015 #include <cmath>
0016
0017
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021
0022 #include "SimG4CMS/HcalTestBeam/interface/HcalTB06Histo.h"
0023
0024
0025
0026 HcalTB06Histo::HcalTB06Histo(const edm::ParameterSet& ps) {
0027 verbose_ = ps.getUntrackedParameter<bool>("Verbose", false);
0028 double em1 = ps.getUntrackedParameter<double>("ETtotMax", 400.);
0029 double em2 = ps.getUntrackedParameter<double>("EHCalMax", 4.0);
0030 mkTree_ = ps.getUntrackedParameter<bool>("MakeTree", false);
0031 eBeam_ = 50.;
0032 mip_ = ps.getParameter<double>("MIP");
0033 edm::LogInfo("HcalTBSim") << "Verbose :" << verbose_ << " MakeTree: " << mkTree_ << " EMax: " << em1 << ":" << em2
0034 << " MIP " << mip_;
0035
0036
0037 edm::Service<TFileService> tfile;
0038
0039 if (!tfile.isAvailable())
0040 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0041 << "please add it to config file";
0042 iniE = tfile->make<TH1D>("iniE", "Incident Energy (GeV)", 4000, 0., em1);
0043 iEta = tfile->make<TH1D>("iEta", "Eta at incidence ", 300, 0., 3.);
0044 iPhi = tfile->make<TH1D>("iPhi", "Phi at incidence ", 300, -1., 1.);
0045 edepS = tfile->make<TH1D>("edepS", "Energy deposit == Total", 4000, 0., em1);
0046 edecS = tfile->make<TH1D>("edecS", "Energy deposit == ECal ", 300, -2., 28.);
0047 edhcS = tfile->make<TH1D>("edhcS", "Energy deposit == HCal ", 4000, 0., em2);
0048 edepN = tfile->make<TH1D>("edepN", "Etot/Ebeam ", 200, -2.5, 2.5);
0049 edecN = tfile->make<TH1D>("edecN", "Eecal/Ebeam ", 200, -2.5, 2.5);
0050 edhcN = tfile->make<TH1D>("edhcN", "Ehcal/Ebeam ", 200, -2.5, 2.5);
0051 emhcN = tfile->make<TH1D>("emhcN", "Ehcal/Ebeam MIP in Ecal", 200, -2.5, 2.5);
0052 edehS = tfile->make<TH2D>("edehS", "Hcal vs Ecal", 100, 0., em1, 100, 0., em2);
0053
0054 if (mkTree_) {
0055 tree_ = tfile->make<TTree>("TB06Sim", "TB06Sim");
0056 tree_->Branch("eBeam_", &eBeam_, "eBeam_/D");
0057 tree_->Branch("etaBeam_", &etaBeam_, "etaBeam_/D");
0058 tree_->Branch("phiBeam_", &phiBeam_, "phiBeam_/D");
0059 tree_->Branch("edepEC_", &edepEC_, "edepEC_/D");
0060 tree_->Branch("edepHB_", &edepHB_, "edepHB_/D");
0061 tree_->Branch("edepHO_", &edepHO_, "edepHO_/D");
0062 tree_->Branch("noiseEC_", &noiseEC_, "noiseEC_/D");
0063 tree_->Branch("noiseHB_", &noiseHB_, "noiseHB_/D");
0064 tree_->Branch("noiseHO_", &noiseHO_, "noiseHO_/D");
0065 tree_->Branch("edepS1_", &edepS1_, "edepS1_/D");
0066 tree_->Branch("edepS2_", &edepS2_, "edepS2_/D");
0067 tree_->Branch("edepS3_", &edepS3_, "edepS3_/D");
0068 tree_->Branch("edepS4_", &edepS4_, "edepS4_/D");
0069 tree_->Branch("edepVC_", &edepVC_, "edepVC_/D");
0070 tree_->Branch("edepS7_", &edepS7_, "edepS7_/D");
0071 tree_->Branch("edepS8_", &edepS8_, "edepS8_/D");
0072 }
0073 }
0074
0075 HcalTB06Histo::~HcalTB06Histo() {}
0076
0077
0078
0079
0080
0081 void HcalTB06Histo::fillPrimary(double energy, double eta, double phi) {
0082 if (verbose_)
0083 edm::LogInfo("HcalTBSim") << "HcalTB06Histo::fillPrimary: Energy " << energy << " Eta " << eta << " Phi " << phi;
0084 eBeam_ = energy;
0085 etaBeam_ = eta;
0086 phiBeam_ = phi;
0087 iniE->Fill(energy);
0088 iEta->Fill(eta);
0089 iPhi->Fill(phi);
0090 }
0091
0092 void HcalTB06Histo::fillEdep(double etots, double eecals, double ehcals) {
0093 if (verbose_)
0094 edm::LogInfo("HcalTBSim") << "HcalTB06Histo:::fillEdep: Simulated Total " << etots << " ECal " << eecals << " HCal "
0095 << ehcals;
0096 edepS->Fill(etots);
0097 edecS->Fill(eecals);
0098 edhcS->Fill(ehcals);
0099 edepN->Fill(etots / eBeam_);
0100 edecN->Fill(eecals / eBeam_);
0101 edhcN->Fill(ehcals / eBeam_);
0102 if (eecals <= mip_) {
0103 emhcN->Fill(etots / eBeam_);
0104 }
0105 edehS->Fill(eecals, ehcals);
0106 }
0107
0108 void HcalTB06Histo::fillTree(std::vector<double>& ecalo, std::vector<double>& etrig) {
0109 if (mkTree_) {
0110 edepEC_ = ecalo[0];
0111 noiseEC_ = ecalo[1];
0112 edepHB_ = ecalo[2];
0113 noiseHB_ = ecalo[3];
0114 edepHO_ = ecalo[4];
0115 noiseHO_ = ecalo[5];
0116 edepS1_ = etrig[0];
0117 edepS2_ = etrig[1];
0118 edepS3_ = etrig[2];
0119 edepS4_ = etrig[3];
0120 edepVC_ = etrig[4];
0121 edepS7_ = etrig[5];
0122 edepS8_ = etrig[6];
0123 tree_->Fill();
0124 if (verbose_)
0125 edm::LogInfo("HcalTBSim") << "HcalTB06Histo:::fillTree: Energies " << edepEC_ << ":" << noiseEC_ << ":" << edepHB_
0126 << ":" << noiseHB_ << ":" << edepHO_ << ":" << noiseHO_ << " Trigger counters "
0127 << edepS1_ << ":" << edepS2_ << ":" << edepS3_ << ":" << edepS4_ << ":" << edepVC_
0128 << ":" << edepS7_ << ":" << edepS8_;
0129 }
0130 }