File indexing completed on 2024-04-06 12:30:12
0001 #include <memory>
0002 #include <string>
0003 #include <fstream>
0004 #include <utility>
0005 #include <vector>
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "SimDataFormats/CaloHit/interface/HFShowerPhoton.h"
0017 #include "SimDataFormats/CaloHit/interface/HFShowerLibraryEventInfo.h"
0018 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020
0021 #include "TFile.h"
0022 #include "TTree.h"
0023
0024 class HcalForwardLibWriter : public edm::one::EDAnalyzer<> {
0025 public:
0026 struct FileHandle {
0027 std::string name;
0028 std::string id;
0029 int momentum;
0030 };
0031 explicit HcalForwardLibWriter(const edm::ParameterSet&);
0032 ~HcalForwardLibWriter() override {}
0033 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0034
0035 private:
0036 void beginJob() override;
0037 void analyze(const edm::Event&, const edm::EventSetup&) override;
0038 void endJob() override;
0039 int readUserData();
0040 int nbins;
0041 int nshowers;
0042 int bsize;
0043 int splitlevel;
0044 int compressionAlgo;
0045 int compressionLevel;
0046
0047 TFile* theFile;
0048 TTree* theTree;
0049 TFile* LibFile;
0050 TTree* LibTree;
0051 TBranch* emBranch;
0052 TBranch* hadBranch;
0053 std::vector<float>* partsEm;
0054 std::vector<float>* partsHad;
0055
0056 std::string theDataFile;
0057 std::vector<FileHandle> theFileHandle;
0058
0059 HFShowerLibraryEventInfo evtInfo;
0060 HFShowerPhotonCollection emColl;
0061 HFShowerPhotonCollection hadColl;
0062 };
0063
0064 HcalForwardLibWriter::HcalForwardLibWriter(const edm::ParameterSet& iConfig) {
0065 edm::ParameterSet theParms = iConfig.getParameter<edm::ParameterSet>("hcalForwardLibWriterParameters");
0066 edm::FileInPath fp = theParms.getParameter<edm::FileInPath>("FileName");
0067 nbins = theParms.getParameter<int>("Nbins");
0068 nshowers = theParms.getParameter<int>("Nshowers");
0069 bsize = theParms.getParameter<int>("BufSize");
0070 splitlevel = theParms.getParameter<int>("SplitLevel");
0071 compressionAlgo = theParms.getParameter<int>("CompressionAlgo");
0072 compressionLevel = theParms.getParameter<int>("CompressionLevel");
0073
0074 std::string pName = fp.fullPath();
0075 if (pName.find('.') == 0)
0076 pName.erase(0, 2);
0077 theDataFile = pName;
0078 readUserData();
0079
0080 edm::Service<TFileService> fs;
0081 fs->file().cd();
0082 fs->file().SetCompressionAlgorithm(compressionAlgo);
0083 fs->file().SetCompressionLevel(compressionLevel);
0084
0085 LibTree = new TTree("HFSimHits", "HFSimHits");
0086
0087
0088
0089 partsEm = new std::vector<float>();
0090 partsHad = new std::vector<float>();
0091 emBranch = LibTree->Branch("emParticles", &partsEm, bsize, splitlevel);
0092 hadBranch = LibTree->Branch("hadParticles", &partsHad, bsize, splitlevel);
0093 }
0094
0095 void HcalForwardLibWriter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0096 edm::ParameterSetDescription desc;
0097 desc.add<edm::FileInPath>("FileName", edm::FileInPath("SimG4CMS/ShowerLibraryProducer/data/fileList.txt"));
0098 desc.add<int>("Nbins", 16);
0099 desc.add<int>("Nshowers", 10000);
0100 desc.add<int>("BufSize", 1);
0101 desc.add<int>("SplitLevel", 2);
0102 desc.add<int>("CompressionAlgo", 4);
0103 desc.add<int>("CompressionLevel", 4);
0104 descriptions.add("hcalForwardLibWriterParameters", desc);
0105 }
0106
0107 void HcalForwardLibWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0108
0109 std::vector<double> en;
0110 double momBin[16] = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
0111 en.reserve(nbins);
0112 for (int i = 0; i < nbins; ++i)
0113 en.push_back(momBin[i]);
0114
0115 int nem = 0;
0116 int nhad = 0;
0117
0118
0119 int n = theFileHandle.size();
0120
0121 for (int i = 0; i < n; ++i) {
0122 std::string fn = theFileHandle[i].name;
0123 std::string particle = theFileHandle[i].id;
0124
0125
0126
0127 TFile* theFile = new TFile(fn.c_str(), "READ");
0128 TTree* theTree = (TTree*)gDirectory->Get("g4SimHits/CherenkovPhotons");
0129 int nphot = 0;
0130
0131 const int size = 10000;
0132 if (nshowers > size) {
0133 edm::LogError("HcalForwardLibWriter") << "Too big Nshowers number";
0134 return;
0135 }
0136
0137 float x[size] = {0.};
0138 float y[size] = {0.};
0139 float z[size] = {0.};
0140 float t[size] = {0.};
0141 float lambda[size] = {0.};
0142 int fiberId[size] = {0};
0143 float primZ;
0144
0145 theTree->SetBranchAddress("nphot", &nphot);
0146 theTree->SetBranchAddress("x", &x);
0147 theTree->SetBranchAddress("y", &y);
0148 theTree->SetBranchAddress("z", &z);
0149 theTree->SetBranchAddress("t", &t);
0150 theTree->SetBranchAddress("lambda", &lambda);
0151 theTree->SetBranchAddress("fiberId", &fiberId);
0152 theTree->SetBranchAddress("primZ", &primZ);
0153 int nentries = int(theTree->GetEntries());
0154 int ngood = 0;
0155
0156 for (int iev = 0; iev < nentries; iev++) {
0157 if (primZ < 990.)
0158 continue;
0159 ngood++;
0160 if (ngood > nshowers)
0161 continue;
0162 unsigned int nph = nphot;
0163 if (particle == "electron") {
0164 emColl.clear();
0165 partsEm->clear();
0166 partsEm->resize(5 * nph);
0167 } else {
0168 hadColl.clear();
0169 partsHad->clear();
0170 partsHad->resize(5 * nph);
0171 }
0172
0173 for (int iph = 0; iph < nphot; ++iph) {
0174 if (fiberId[iph] != 1) {
0175 z[iph] = -z[iph];
0176 }
0177
0178 if (particle == "electron") {
0179 (*partsEm)[iph] = (x[iph]);
0180 (*partsEm)[iph + 1 * nph] = (y[iph]);
0181 (*partsEm)[iph + 2 * nph] = (z[iph]);
0182 (*partsEm)[iph + 3 * nph] = (t[iph]);
0183 (*partsEm)[iph + 4 * nph] = (lambda[iph]);
0184 } else {
0185 (*partsHad)[iph] = (x[iph]);
0186 (*partsHad)[iph + 1 * nph] = (y[iph]);
0187 (*partsHad)[iph + 2 * nph] = (z[iph]);
0188 (*partsHad)[iph + 3 * nph] = (t[iph]);
0189 (*partsHad)[iph + 4 * nph] = (lambda[iph]);
0190 }
0191 }
0192
0193
0194 if (particle == "electron") {
0195 LibTree->SetEntries(nem);
0196 emBranch->Fill();
0197 nem++;
0198 emColl.clear();
0199 } else {
0200 LibTree->SetEntries(nhad);
0201 nhad++;
0202 hadBranch->Fill();
0203 hadColl.clear();
0204 }
0205 }
0206
0207 theFile->Close();
0208 }
0209
0210 }
0211
0212 void HcalForwardLibWriter::beginJob() {}
0213
0214 void HcalForwardLibWriter::endJob() {
0215 edm::Service<TFileService> fs;
0216 fs->file().cd();
0217 LibTree->Write();
0218 LibTree->Print();
0219 }
0220
0221 int HcalForwardLibWriter::readUserData(void) {
0222 std::ifstream input(theDataFile.c_str());
0223 if (input.fail()) {
0224 return 0;
0225 }
0226 std::string theFileName, thePID;
0227 int mom;
0228 int k = 0;
0229 while (!input.eof()) {
0230 input >> theFileName >> thePID >> mom;
0231 if (!input.fail()) {
0232 FileHandle aFile;
0233 aFile.name = theFileName;
0234 aFile.id = thePID;
0235 aFile.momentum = mom;
0236 theFileHandle.push_back(aFile);
0237 ++k;
0238 } else {
0239 input.clear();
0240 }
0241 input.ignore(999, '\n');
0242 }
0243 return k;
0244 }
0245
0246
0247 DEFINE_FWK_MODULE(HcalForwardLibWriter);