File indexing completed on 2024-05-10 02:21:22
0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "FWCore/Utilities/interface/InputTag.h"
0012
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0015
0016 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
0017 #include "SimDataFormats/CaloHit/interface/HFShowerPhoton.h"
0018
0019 #include <CLHEP/Units/SystemOfUnits.h>
0020 #include "TFile.h"
0021 #include "TTree.h"
0022 #include "TBranch.h"
0023 #include "TH1F.h"
0024
0025 #include <string>
0026 #include <vector>
0027
0028 class HFShowerLibraryAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0029 public:
0030 HFShowerLibraryAnalyzer(const edm::ParameterSet& ps);
0031 ~HFShowerLibraryAnalyzer() override;
0032 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0033
0034 protected:
0035 void beginJob() override {}
0036 void analyze(edm::Event const&, edm::EventSetup const&) override {}
0037 void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0038 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0039
0040 private:
0041 void bookHistos();
0042 void getRecord(int, int);
0043 void loadEventInfo(TBranch*);
0044
0045 TFile* hf_;
0046 TBranch *emBranch_, *hadBranch_;
0047 bool verbose_, newForm_, v3version_;
0048 int nMomBin_, totEvents_, evtPerBin_;
0049 float libVers_, listVersion_;
0050 std::vector<double> pmom_;
0051 HFShowerPhotonCollection* photo_;
0052 HFShowerPhotonCollection photon_;
0053 std::vector<TH1F*> h_x_[2], h_y_[2], h_z_[2], h_t_[2], h_l_[2], h_nl_[2], h_ns_[2];
0054 };
0055
0056 HFShowerLibraryAnalyzer::HFShowerLibraryAnalyzer(edm::ParameterSet const& ps)
0057 : hf_(nullptr), emBranch_(nullptr), hadBranch_(nullptr) {
0058 usesResource(TFileService::kSharedResource);
0059
0060 edm::FileInPath fp("SimG4CMS/Calo/data/" + ps.getParameter<std::string>("FileName"));
0061 std::string pTreeName = fp.fullPath();
0062 std::string emName = ps.getParameter<std::string>("TreeEMID");
0063 std::string hadName = ps.getParameter<std::string>("TreeHadID");
0064 std::string branchEvInfo = ps.getParameter<std::string>("BranchEvt");
0065 std::string branchPre = ps.getParameter<std::string>("BranchPre");
0066 std::string branchPost = ps.getParameter<std::string>("BranchPost");
0067 verbose_ = ps.getParameter<bool>("Verbosity");
0068 evtPerBin_ = ps.getParameter<bool>("EventPerBin");
0069
0070 if (pTreeName.find('.') == 0)
0071 pTreeName.erase(0, 2);
0072 const char* nTree = pTreeName.c_str();
0073 hf_ = TFile::Open(nTree);
0074
0075 if (hf_->IsOpen()) {
0076 edm::LogVerbatim("HFShower") << "HFShowerLibrary: opening " << nTree << " successfully";
0077 newForm_ = (branchEvInfo.empty());
0078 TTree* event = ((newForm_) ? static_cast<TTree*>(hf_->Get("HFSimHits")) : static_cast<TTree*>(hf_->Get("Events")));
0079 if (event) {
0080 TBranch* evtInfo(nullptr);
0081 if (!newForm_) {
0082 std::string info = branchEvInfo + branchPost;
0083 evtInfo = event->GetBranch(info.c_str());
0084 }
0085 if (evtInfo || newForm_) {
0086 loadEventInfo(evtInfo);
0087 } else {
0088 edm::LogError("HFShower") << "HFShowerLibrayEventInfo Branch does not exist in Event";
0089 throw cms::Exception("Unknown", "HFShowerLibraryAnalyzer") << "Event information absent\n";
0090 }
0091 } else {
0092 edm::LogError("HFShower") << "Events Tree does not exist";
0093 throw cms::Exception("Unknown", "HFShowerLibraryAnalyzer") << "Events tree absent\n";
0094 }
0095
0096 std::stringstream ss;
0097 ss << "HFShowerLibraryAnalyzer: Library " << libVers_ << " ListVersion " << listVersion_ << " Events Total "
0098 << totEvents_ << " and " << evtPerBin_ << " per bin\n";
0099 ss << "HFShowerLibraryAnalyzer: Energies (GeV) with " << nMomBin_ << " bins\n";
0100 for (int i = 0; i < nMomBin_; ++i) {
0101 if (i / 10 * 10 == i && i > 0) {
0102 ss << "\n";
0103 }
0104 ss << " " << pmom_[i] / CLHEP::GeV;
0105 }
0106 edm::LogVerbatim("HFShower") << ss.str();
0107
0108 std::string nameBr = branchPre + emName + branchPost;
0109 emBranch_ = event->GetBranch(nameBr.c_str());
0110 if (verbose_)
0111 emBranch_->Print();
0112 nameBr = branchPre + hadName + branchPost;
0113 hadBranch_ = event->GetBranch(nameBr.c_str());
0114 if (verbose_)
0115 hadBranch_->Print();
0116
0117 v3version_ = (emBranch_->GetClassName() == std::string("vector<float>")) ? true : false;
0118 edm::LogVerbatim("HFShower")
0119 << " HFShowerLibraryAnalyzer:Branch " << emName << " has " << emBranch_->GetEntries() << " entries and Branch "
0120 << hadName << " has " << hadBranch_->GetEntries()
0121 << " entries\n HFShowerLibraryAnalyzer::No packing information - Assume x, y, z are not in packed form";
0122 photo_ = new HFShowerPhotonCollection;
0123
0124 bookHistos();
0125 } else {
0126 edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree << " failed";
0127 }
0128 }
0129
0130 HFShowerLibraryAnalyzer::~HFShowerLibraryAnalyzer() {
0131 if (hf_)
0132 hf_->Close();
0133 delete photo_;
0134 }
0135
0136 void HFShowerLibraryAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0137 edm::ParameterSetDescription desc;
0138 desc.add<std::string>("FileName", "HFShowerLibrary_10000.root");
0139 desc.add<std::string>("TreeEMID", "emParticles");
0140 desc.add<std::string>("TreeHadID", "hadParticles");
0141 desc.add<std::string>("BranchEvt", "");
0142 desc.add<std::string>("BranchPre", "");
0143 desc.add<std::string>("BranchPost", "");
0144 desc.add<bool>("Verbosity", false);
0145 desc.add<int>("EventPerBin", 10000);
0146 descriptions.add("hfShowerLibaryAnalysis", desc);
0147 }
0148
0149 void HFShowerLibraryAnalyzer::bookHistos() {
0150 edm::Service<TFileService> tfile;
0151 if (!tfile.isAvailable())
0152 throw cms::Exception("Unknown", "HFShowerLibraryAnalyzer")
0153 << "TFileService unavailable: please add it to config file";
0154 char name[20], title[40], titlx[120];
0155 TH1F* hist;
0156 for (int i = 0; i < 2; ++i) {
0157 std::string type = (i == 0) ? "EM" : "Hadron";
0158 for (int k = 0; k < nMomBin_; ++k) {
0159 sprintf(title, "Showers for p = %6.1f GeV", pmom_[k] / CLHEP::GeV);
0160 sprintf(name, "X%d%d", k, i);
0161 sprintf(titlx, "x coordinate of %s shower (mm)", type.c_str());
0162 hist = tfile->make<TH1F>(name, title, 200, -500.0, 500.0);
0163 hist->GetXaxis()->SetTitle(titlx);
0164 hist->GetYaxis()->SetTitle("Shower");
0165 hist->GetYaxis()->SetTitleOffset(1.2);
0166 hist->Sumw2();
0167 h_x_[i].emplace_back(hist);
0168 sprintf(name, "Y%d%d", k, i);
0169 sprintf(titlx, "y coordinate of %s shower (mm)", type.c_str());
0170 hist = tfile->make<TH1F>(name, title, 200, -500.0, 500.0);
0171 hist->GetXaxis()->SetTitle(titlx);
0172 hist->GetYaxis()->SetTitle("Shower");
0173 hist->GetYaxis()->SetTitleOffset(1.2);
0174 hist->Sumw2();
0175 h_y_[i].emplace_back(hist);
0176 sprintf(name, "Z%d%d", k, i);
0177 sprintf(titlx, "z coordinate of %s shower (mm)", type.c_str());
0178 hist = tfile->make<TH1F>(name, title, 200, -2000.0, 2000.0);
0179 hist->GetXaxis()->SetTitle(titlx);
0180 hist->GetYaxis()->SetTitle("Shower");
0181 hist->GetYaxis()->SetTitleOffset(1.2);
0182 hist->Sumw2();
0183 h_z_[i].emplace_back(hist);
0184 sprintf(name, "T%d%d", k, i);
0185 sprintf(titlx, "Time of %s shower (ns)", type.c_str());
0186 hist = tfile->make<TH1F>(name, title, 200, 0.0, 50.0);
0187 hist->GetXaxis()->SetTitle(titlx);
0188 hist->GetYaxis()->SetTitle("Shower");
0189 hist->GetYaxis()->SetTitleOffset(1.2);
0190 hist->Sumw2();
0191 h_t_[i].emplace_back(hist);
0192 sprintf(name, "L%d%d", k, i);
0193 sprintf(titlx, "Lambda of %s shower photon (nm)", type.c_str());
0194 hist = tfile->make<TH1F>(name, title, 200, 300.0, 800.0);
0195 hist->GetXaxis()->SetTitle(titlx);
0196 hist->GetYaxis()->SetTitle("Shower");
0197 hist->GetYaxis()->SetTitleOffset(1.2);
0198 h_l_[i].emplace_back(hist);
0199 hist->Sumw2();
0200 sprintf(name, "NL%d%d", k, i);
0201 sprintf(titlx, "Number of %s shower photons in long fibers", type.c_str());
0202 hist = tfile->make<TH1F>(name, title, 3000, 0.0, 3000.0);
0203 hist->GetXaxis()->SetTitle(titlx);
0204 hist->GetYaxis()->SetTitle("Shower");
0205 hist->GetYaxis()->SetTitleOffset(1.2);
0206 h_nl_[i].emplace_back(hist);
0207 sprintf(name, "NS%d%d", k, i);
0208 sprintf(titlx, "Number of %s shower photons in short fibes", type.c_str());
0209 hist = tfile->make<TH1F>(name, title, 3000, 0.0, 3000.0);
0210 hist->GetXaxis()->SetTitle(titlx);
0211 hist->GetYaxis()->SetTitle("Shower");
0212 hist->GetYaxis()->SetTitleOffset(1.2);
0213 h_ns_[i].emplace_back(hist);
0214 }
0215 }
0216
0217
0218 for (int type = 0; type < 2; ++type) {
0219 for (int k = 0; k < nMomBin_; ++k) {
0220 for (int j = 0; j < evtPerBin_; ++j) {
0221 int irc = k * evtPerBin_ + j + 1;
0222 getRecord(type, irc);
0223 int nPhoton = (newForm_) ? photo_->size() : photon_.size();
0224 int nlong = 0, nshort = 0;
0225 for (int i = 0; i < nPhoton; i++) {
0226 if (newForm_) {
0227 if (photo_->at(i).z() > 0) {
0228 nlong++;
0229 } else {
0230 nshort++;
0231 }
0232 h_x_[type][k]->Fill((photo_->at(i)).x());
0233 h_y_[type][k]->Fill((photo_->at(i)).y());
0234 h_z_[type][k]->Fill((photo_->at(i)).z());
0235 h_t_[type][k]->Fill((photo_->at(i)).t());
0236 h_l_[type][k]->Fill((photo_->at(i)).lambda());
0237 } else {
0238 if (photon_[i].z() > 0) {
0239 nlong++;
0240 } else {
0241 nshort++;
0242 }
0243 h_x_[type][k]->Fill(photon_[i].x());
0244 h_y_[type][k]->Fill(photon_[i].y());
0245 h_z_[type][k]->Fill(photon_[i].z());
0246 h_t_[type][k]->Fill(photon_[i].t());
0247 h_l_[type][k]->Fill(photon_[i].lambda());
0248 }
0249 }
0250 h_nl_[type][k]->Fill(nlong);
0251 h_ns_[type][k]->Fill(nshort);
0252 }
0253 }
0254 }
0255 }
0256
0257 void HFShowerLibraryAnalyzer::getRecord(int type, int record) {
0258 int nrc = record - 1;
0259 photon_.clear();
0260 photo_->clear();
0261 if (type > 0) {
0262 if (newForm_) {
0263 if (!v3version_) {
0264 hadBranch_->SetAddress(&photo_);
0265 hadBranch_->GetEntry(nrc + totEvents_);
0266 } else {
0267 std::vector<float> t;
0268 std::vector<float>* tp = &t;
0269 hadBranch_->SetAddress(&tp);
0270 hadBranch_->GetEntry(nrc + totEvents_);
0271 unsigned int tSize = t.size() / 5;
0272 photo_->reserve(tSize);
0273 for (unsigned int i = 0; i < tSize; i++) {
0274 photo_->push_back(
0275 HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
0276 }
0277 }
0278 } else {
0279 hadBranch_->SetAddress(&photon_);
0280 hadBranch_->GetEntry(nrc);
0281 }
0282 } else {
0283 if (newForm_) {
0284 if (!v3version_) {
0285 emBranch_->SetAddress(&photo_);
0286 emBranch_->GetEntry(nrc);
0287 } else {
0288 std::vector<float> t;
0289 std::vector<float>* tp = &t;
0290 emBranch_->SetAddress(&tp);
0291 emBranch_->GetEntry(nrc);
0292 unsigned int tSize = t.size() / 5;
0293 photo_->reserve(tSize);
0294 for (unsigned int i = 0; i < tSize; i++) {
0295 photo_->push_back(
0296 HFShowerPhoton(t[i], t[1 * tSize + i], t[2 * tSize + i], t[3 * tSize + i], t[4 * tSize + i]));
0297 }
0298 }
0299 } else {
0300 emBranch_->SetAddress(&photon_);
0301 emBranch_->GetEntry(nrc);
0302 }
0303 }
0304 if (verbose_) {
0305 int nPhoton = (newForm_) ? photo_->size() : photon_.size();
0306 edm::LogVerbatim("HFShower") << "getRecord: Record " << record << " of type " << type << " with " << nPhoton
0307 << " photons";
0308 for (int j = 0; j < nPhoton; j++)
0309 if (newForm_)
0310 edm::LogVerbatim("HFShower") << "Photon " << j << " " << photo_->at(j);
0311 else
0312 edm::LogVerbatim("HFShower") << "Photon " << j << " " << photon_[j];
0313 }
0314 }
0315
0316 void HFShowerLibraryAnalyzer::loadEventInfo(TBranch* branch) {
0317 if (branch) {
0318 std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
0319 branch->SetAddress(&eventInfoCollection);
0320 branch->GetEntry(0);
0321 edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo Collection of size "
0322 << eventInfoCollection.size() << " records";
0323 totEvents_ = eventInfoCollection[0].totalEvents();
0324 nMomBin_ = eventInfoCollection[0].numberOfBins();
0325 evtPerBin_ = eventInfoCollection[0].eventsPerBin();
0326 libVers_ = eventInfoCollection[0].showerLibraryVersion();
0327 listVersion_ = eventInfoCollection[0].physListVersion();
0328 pmom_ = eventInfoCollection[0].energyBins();
0329 } else {
0330 edm::LogVerbatim("HFShower") << "HFShowerLibrary::loadEventInfo loads EventInfo from hardwired numbers";
0331 nMomBin_ = 16;
0332 totEvents_ = nMomBin_ * evtPerBin_;
0333 libVers_ = 1.1;
0334 listVersion_ = 3.6;
0335 pmom_ = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
0336 }
0337 for (int i = 0; i < nMomBin_; i++)
0338 pmom_[i] *= CLHEP::GeV;
0339 }
0340
0341
0342 DEFINE_FWK_MODULE(HFShowerLibraryAnalyzer);
