File indexing completed on 2024-04-06 12:30:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <vector>
0022 #include <string>
0023
0024 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0025
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0033 #include "FWCore/ServiceRegistry/interface/Service.h"
0034 #include "FWCore/Utilities/interface/InputTag.h"
0035
0036 #include "SimDataFormats/CaloHit/interface/HFShowerPhoton.h"
0037 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
0038
0039 #include "TFile.h"
0040 #include "TTree.h"
0041 #include "TBranch.h"
0042 #include "TH1F.h"
0043 #include "TH1I.h"
0044
0045 class AnalyzeTuples : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0046 public:
0047 explicit AnalyzeTuples(const edm::ParameterSet&);
0048 ~AnalyzeTuples() override = default;
0049 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0050
0051 private:
0052 void beginJob() override;
0053 void analyze(const edm::Event&, const edm::EventSetup&) override;
0054 void endJob() override;
0055 void loadEventInfo(TBranch*);
0056 void getRecord(int type, int record);
0057 TFile* hf;
0058 TBranch *emBranch, *hadBranch;
0059
0060 int nMomBin, totEvents, evtPerBin;
0061 float libVers, listVersion;
0062 std::vector<double> pmom;
0063 std::vector<HFShowerPhoton> photon;
0064
0065 TH1I* hNPELongElec[12];
0066 TH1I* hNPEShortElec[12];
0067 TH1I* hNPELongPion[12];
0068 TH1I* hNPEShortPion[12];
0069 };
0070
0071 AnalyzeTuples::AnalyzeTuples(const edm::ParameterSet& iConfig) {
0072 usesResource(TFileService::kSharedResource);
0073
0074 edm::LogVerbatim("HcalSim") << "analyzetuples a buraya girdi";
0075 edm::FileInPath fp = iConfig.getParameter<edm::FileInPath>("FileName");
0076 std::string pTreeName = fp.fullPath();
0077 std::string emName = iConfig.getParameter<std::string>("TreeEMID");
0078 std::string hadName = iConfig.getParameter<std::string>("TreeHadID");
0079 std::string branchEvInfo = iConfig.getUntrackedParameter<std::string>(
0080 "BranchEvt", "HFShowerLibraryEventInfos_hfshowerlib_HFShowerLibraryEventInfo");
0081 std::string branchPre = iConfig.getUntrackedParameter<std::string>("BranchPre", "HFShowerPhotons_hfshowerlib_");
0082 std::string branchPost = iConfig.getUntrackedParameter<std::string>("BranchPost", "_R.obj");
0083
0084 if (pTreeName.find(".") == 0)
0085 pTreeName.erase(0, 2);
0086 const char* nTree = pTreeName.c_str();
0087 hf = TFile::Open(nTree);
0088
0089 if (!hf->IsOpen()) {
0090 edm::LogError("HFShower") << "HFShowerLibrary: opening " << nTree << " failed";
0091 throw cms::Exception("Unknown", "HFShowerLibrary") << "Opening of " << pTreeName << " fails\n";
0092 } else {
0093 edm::LogInfo("HFShower") << "HFShowerLibrary: opening " << nTree << " successfully";
0094 }
0095
0096 TTree* event = (TTree*)hf->Get("Events");
0097 if (event) {
0098 std::string info = branchEvInfo + branchPost;
0099 TBranch* evtInfo = event->GetBranch(info.c_str());
0100 if (evtInfo) {
0101 loadEventInfo(evtInfo);
0102 } else {
0103 edm::LogError("HFShower") << "HFShowerLibrary: HFShowerLibrayEventInfo"
0104 << " Branch does not exist in Event";
0105 throw cms::Exception("Unknown", "HFShowerLibrary") << "Event information absent\n";
0106 }
0107 } else {
0108 edm::LogError("HFShower") << "HFShowerLibrary: Events Tree does not "
0109 << "exist";
0110 throw cms::Exception("Unknown", "HFShowerLibrary") << "Events tree absent\n";
0111 }
0112
0113 edm::LogInfo("HFShower") << "HFShowerLibrary: Library " << libVers << " ListVersion " << listVersion
0114 << " Events Total " << totEvents << " and " << evtPerBin << " per bin";
0115 edm::LogInfo("HFShower") << "HFShowerLibrary: Energies (GeV) with " << nMomBin << " bins";
0116 for (int i = 0; i < nMomBin; i++)
0117 edm::LogInfo("HFShower") << "HFShowerLibrary: pmom[" << i << "] = " << pmom[i] << " GeV";
0118
0119 std::string nameBr = branchPre + emName + branchPost;
0120 emBranch = event->GetBranch(nameBr.c_str());
0121 nameBr = branchPre + hadName + branchPost;
0122 hadBranch = event->GetBranch(nameBr.c_str());
0123 edm::LogInfo("HFShower") << "HFShowerLibrary:Branch " << emName << " has " << emBranch->GetEntries()
0124 << " entries and Branch " << hadName << " has " << hadBranch->GetEntries() << " entries";
0125 edm::LogInfo("HFShower") << "HFShowerLibrary::No packing information -"
0126 << " Assume x, y, z are not in packed form";
0127 }
0128
0129 void AnalyzeTuples::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0130 edm::ParameterSetDescription desc;
0131 desc.add<std::string>("FileName", "SimG4CMS/Calo/data/HFShowerLibrary_oldpmt_noatt_eta4_16en_v3.root");
0132 desc.add<std::string>("TreeEMID", "emParticles");
0133 desc.add<std::string>("TreeHadID", "hadParticles");
0134 desc.add<std::string>("BranchEvt", "");
0135 desc.add<std::string>("BranchPre", "");
0136 desc.add<std::string>("BranchPost", "");
0137 descriptions.add("analyzeTuples", desc);
0138 }
0139
0140 void AnalyzeTuples::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0141 for (int ibin = 0; ibin < 12; ++ibin) {
0142 int min = evtPerBin * (ibin);
0143 int max = evtPerBin * (ibin + 1);
0144 for (int i = min; i < max; ++i) {
0145 getRecord(0, i);
0146 int npe_long = 0;
0147 int npe_short = 0;
0148 edm::LogVerbatim("HcalSim") << "phptons size" << photon.size();
0149 for (int j = 0; j < int(photon.size()); ++j) {
0150
0151 if (photon[j].z() < 0) {
0152
0153 ++npe_short;
0154 } else {
0155
0156 ++npe_long;
0157 edm::LogVerbatim("HcalSim") << photon[j].z();
0158 }
0159 }
0160 hNPELongElec[ibin]->Fill(npe_long);
0161 edm::LogVerbatim("HcalSim") << ibin << npe_long;
0162 hNPEShortElec[ibin]->Fill(npe_short);
0163 }
0164 }
0165 for (int ibin = 0; ibin < 12; ++ibin) {
0166 int min = evtPerBin * (ibin);
0167 int max = evtPerBin * (ibin + 1);
0168 for (int i = min; i < max; ++i) {
0169 getRecord(1, i);
0170 int npe_long = 0;
0171 int npe_short = 0;
0172 for (int j = 0; j < int(photon.size()); ++j) {
0173
0174 if (photon[j].z() < 0) {
0175
0176 ++npe_short;
0177 } else {
0178
0179 ++npe_long;
0180 }
0181 }
0182 hNPELongPion[ibin]->Fill(npe_long);
0183 hNPEShortPion[ibin]->Fill(npe_short);
0184 }
0185 }
0186 }
0187
0188 void AnalyzeTuples::beginJob() {
0189 edm::Service<TFileService> fs;
0190 TFileDirectory HFDir = fs->mkdir("HF");
0191 char title[128];
0192 for (int i = 0; i < int(pmom.size()); ++i) {
0193 sprintf(title, "NPELongElec_Mom_%i", int(pmom[i]));
0194 int maxBin = int(pmom[i] + 50);
0195 hNPELongElec[i] = HFDir.make<TH1I>(title, "NPE Long", 140, 0, 140);
0196 sprintf(title, "NPEShortElec_Mom_%i", int(pmom[i]));
0197 hNPEShortElec[i] = HFDir.make<TH1I>(title, "NPE Short", maxBin, 0, maxBin);
0198 sprintf(title, "NPELongPion_Mom_%i", int(pmom[i]));
0199 hNPELongPion[i] = HFDir.make<TH1I>(title, "NPE Long", maxBin, 0, maxBin);
0200 sprintf(title, "NPEShortPion_Mom_%i", int(pmom[i]));
0201 hNPEShortPion[i] = HFDir.make<TH1I>(title, "NPE Short", maxBin, 0, maxBin);
0202 }
0203 }
0204
0205 void AnalyzeTuples::endJob() {}
0206
0207 void AnalyzeTuples::loadEventInfo(TBranch* branch) {
0208 std::vector<HFShowerLibraryEventInfo> eventInfoCollection;
0209 branch->SetAddress(&eventInfoCollection);
0210 branch->GetEntry(0);
0211 edm::LogInfo("HFShower") << "HFShowerLibrary::loadEventInfo loads "
0212 << " EventInfo Collection of size " << eventInfoCollection.size() << " records";
0213 totEvents = eventInfoCollection[0].totalEvents();
0214 nMomBin = eventInfoCollection[0].numberOfBins();
0215 evtPerBin = eventInfoCollection[0].eventsPerBin();
0216 libVers = eventInfoCollection[0].showerLibraryVersion();
0217 listVersion = eventInfoCollection[0].physListVersion();
0218 pmom = eventInfoCollection[0].energyBins();
0219 }
0220
0221 void AnalyzeTuples::getRecord(int type, int record) {
0222 int nrc = record - 1;
0223 photon.clear();
0224 if (type > 0) {
0225 hadBranch->SetAddress(&photon);
0226 hadBranch->GetEntry(nrc);
0227 } else {
0228 emBranch->SetAddress(&photon);
0229 emBranch->GetEntry(nrc);
0230 }
0231 #ifdef EDM_ML_DEBUG
0232 int nPhoton = photon.size();
0233 edm::LogVerbatim("HFShower") << "HFShowerLibrary::getRecord: Record " << record << " of type " << type << " with "
0234 << nPhoton << " photons";
0235 for (int j = 0; j < nPhoton; j++)
0236 edm::LogVerbatim("HFShower") << "Photon " << j << " " << photon[j];
0237 #endif
0238 }
0239
0240
0241 DEFINE_FWK_MODULE(AnalyzeTuples);
