ShowerLibraryProducer/plugins/HcalForwardLibWriter.cc

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   //https://root.cern/root/html534/TTree.html
0088   // TBranch*Branch(const char* name, const char* classname, void** obj, Int_t bufsize = 32000, Int_t splitlevel = 99)
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   // Event info
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   //shower photons
0119   int n = theFileHandle.size();
0120   // cycle over files ++++++++++++++++++++++++++++++++++++++++++++++++++++
0121   for (int i = 0; i < n; ++i) {
0122     std::string fn = theFileHandle[i].name;
0123     std::string particle = theFileHandle[i].id;
0124
0125     //    edm::LogVerbatim("HcalSim") << "*** Input file  " << i << "   " << fn;
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;  // added
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);  // added
0153     int nentries = int(theTree->GetEntries());
0154     int ngood = 0;
0155     // cycle over showers ====================================================
0156     for (int iev = 0; iev < nentries; iev++) {
0157       if (primZ < 990.)
0158         continue;  // exclude showers with interactions in front of HF (1m of air)
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       // cycle over photons in shower -------------------------------------------
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       // end of cycle over photons in shower -------------------------------------------
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     // end of cycle over showers ====================================================
0207     theFile->Close();
0208   }
0209   // end of cycle over files ++++++++++++++++++++++++++++++++++++++++++++++++++++
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 //define this as a plug-in
0247 DEFINE_FWK_MODULE(HcalForwardLibWriter);