File indexing completed on 2023-01-16 23:37:23
0001 #include "SimG4CMS/ShowerLibraryProducer/interface/HcalForwardLibWriter.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 HcalForwardLibWriter::HcalForwardLibWriter(const edm::ParameterSet& iConfig) {
0005 edm::ParameterSet theParms = iConfig.getParameter<edm::ParameterSet>("hcalForwardLibWriterParameters");
0006 edm::FileInPath fp = theParms.getParameter<edm::FileInPath>("FileName");
0007 nbins = theParms.getParameter<int>("Nbins");
0008 nshowers = theParms.getParameter<int>("Nshowers");
0009 bsize = theParms.getParameter<int>("BufSize");
0010 splitlevel = theParms.getParameter<int>("SplitLevel");
0011 compressionAlgo = theParms.getParameter<int>("CompressionAlgo");
0012 compressionLevel = theParms.getParameter<int>("CompressionLevel");
0013
0014 std::string pName = fp.fullPath();
0015 if (pName.find('.') == 0)
0016 pName.erase(0, 2);
0017 theDataFile = pName;
0018 readUserData();
0019
0020 edm::Service<TFileService> fs;
0021 fs->file().cd();
0022 fs->file().SetCompressionAlgorithm(compressionAlgo);
0023 fs->file().SetCompressionLevel(compressionLevel);
0024
0025 LibTree = new TTree("HFSimHits", "HFSimHits");
0026
0027
0028
0029 emBranch = LibTree->Branch("emParticles", "HFShowerPhotons-emParticles", &emColl, bsize, splitlevel);
0030 hadBranch = LibTree->Branch("hadParticles", "HFShowerPhotons-hadParticles", &hadColl, bsize, splitlevel);
0031 }
0032
0033 void HcalForwardLibWriter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0034 edm::ParameterSetDescription desc;
0035 desc.add<edm::FileInPath>("FileName", edm::FileInPath("SimG4CMS/ShowerLibraryProducer/data/fileList.txt"));
0036 desc.add<int>("Nbins", 16);
0037 desc.add<int>("Nshowers", 10000);
0038 desc.add<int>("BufSize", 1);
0039 desc.add<int>("SplitLevel", 2);
0040 desc.add<int>("CompressionAlgo", 4);
0041 desc.add<int>("CompressionLevel", 4);
0042 descriptions.add("hcalForwardLibWriterParameters", desc);
0043 }
0044
0045 void HcalForwardLibWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0046
0047 std::vector<double> en;
0048 double momBin[16] = {2, 3, 5, 7, 10, 15, 20, 30, 50, 75, 100, 150, 250, 350, 500, 1000};
0049 en.reserve(nbins);
0050 for (int i = 0; i < nbins; ++i)
0051 en.push_back(momBin[i]);
0052
0053 int nem = 0;
0054 int nhad = 0;
0055
0056
0057 int n = theFileHandle.size();
0058
0059 for (int i = 0; i < n; ++i) {
0060 std::string fn = theFileHandle[i].name;
0061 std::string particle = theFileHandle[i].id;
0062
0063
0064
0065 TFile* theFile = new TFile(fn.c_str(), "READ");
0066 TTree* theTree = (TTree*)gDirectory->Get("g4SimHits/CherenkovPhotons");
0067 int nphot = 0;
0068
0069 const int size = 10000;
0070 if (nshowers > size) {
0071 edm::LogError("HcalForwardLibWriter") << "Too big Nshowers number";
0072 return;
0073 }
0074
0075 float x[size] = {0.};
0076 float y[size] = {0.};
0077 float z[size] = {0.};
0078 float t[size] = {0.};
0079 float lambda[size] = {0.};
0080 int fiberId[size] = {0};
0081 float primZ;
0082
0083 theTree->SetBranchAddress("nphot", &nphot);
0084 theTree->SetBranchAddress("x", &x);
0085 theTree->SetBranchAddress("y", &y);
0086 theTree->SetBranchAddress("z", &z);
0087 theTree->SetBranchAddress("t", &t);
0088 theTree->SetBranchAddress("lambda", &lambda);
0089 theTree->SetBranchAddress("fiberId", &fiberId);
0090 theTree->SetBranchAddress("primZ", &primZ);
0091 int nentries = int(theTree->GetEntries());
0092 int ngood = 0;
0093
0094 for (int iev = 0; iev < nentries; iev++) {
0095 if (primZ < 990.)
0096 continue;
0097 ngood++;
0098 if (ngood > nshowers)
0099 continue;
0100 if (particle == "electron") {
0101 emColl.clear();
0102 } else {
0103 hadColl.clear();
0104 }
0105 float nphot_long = 0;
0106 float nphot_short = 0;
0107
0108 for (int iph = 0; iph < nphot; ++iph) {
0109 if (fiberId[iph] == 1) {
0110 nphot_long++;
0111 } else {
0112 nphot_short++;
0113 z[iph] = -z[iph];
0114 }
0115
0116 HFShowerPhoton::Point pos(x[iph], y[iph], z[iph]);
0117 HFShowerPhoton aPhoton(pos, t[iph], lambda[iph]);
0118 if (particle == "electron") {
0119 emColl.push_back(aPhoton);
0120 } else {
0121 hadColl.push_back(aPhoton);
0122 }
0123 }
0124
0125
0126 if (particle == "electron") {
0127 LibTree->SetEntries(nem + 1);
0128 emBranch->Fill();
0129 nem++;
0130 emColl.clear();
0131 } else {
0132 LibTree->SetEntries(nhad + 1);
0133 nhad++;
0134 hadBranch->Fill();
0135 hadColl.clear();
0136 }
0137 }
0138
0139 theFile->Close();
0140 }
0141
0142 }
0143
0144 void HcalForwardLibWriter::beginJob() {}
0145
0146 void HcalForwardLibWriter::endJob() {
0147 edm::Service<TFileService> fs;
0148 fs->file().cd();
0149 LibTree->Write();
0150 LibTree->Print();
0151 }
0152
0153 int HcalForwardLibWriter::readUserData(void) {
0154 std::ifstream input(theDataFile.c_str());
0155 if (input.fail()) {
0156 return 0;
0157 }
0158 std::string theFileName, thePID;
0159 int mom;
0160 int k = 0;
0161 while (!input.eof()) {
0162 input >> theFileName >> thePID >> mom;
0163 if (!input.fail()) {
0164 FileHandle aFile;
0165 aFile.name = theFileName;
0166 aFile.id = thePID;
0167 aFile.momentum = mom;
0168 theFileHandle.push_back(aFile);
0169 ++k;
0170 } else {
0171 input.clear();
0172 }
0173 input.ignore(999, '\n');
0174 }
0175 return k;
0176 }
0177
0178
0179 DEFINE_FWK_MODULE(HcalForwardLibWriter);