ShowerLibraryProducer/plugins/HcalForwardLibWriter.cc

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   //https://root.cern/root/html534/TTree.html
0028   // TBranch*Branch(const char* name, const char* classname, void** obj, Int_t bufsize = 32000, Int_t splitlevel = 99)
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   // Event info
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   //shower photons
0057   int n = theFileHandle.size();
0058   // cycle over files ++++++++++++++++++++++++++++++++++++++++++++++++++++
0059   for (int i = 0; i < n; ++i) {
0060     std::string fn = theFileHandle[i].name;
0061     std::string particle = theFileHandle[i].id;
0062
0063     //    edm::LogVerbatim("HcalSim") << "*** Input file  " << i << "   " << fn;
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;  // added
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);  // added
0091     int nentries = int(theTree->GetEntries());
0092     int ngood = 0;
0093     int nbytes = 0;
0094     // cycle over showers ====================================================
0095     for (int iev = 0; iev < nentries; iev++) {
0096       nbytes += theTree->GetEntry(iev);
0097       if (primZ < 990.)
0098         continue;  // exclude showers with interactions in front of HF (1m of air)
0099       ngood++;
0100       if (ngood > nshowers)
0101         continue;
0102       if (particle == "electron") {
0103         emColl.clear();
0104       } else {
0105         hadColl.clear();
0106       }
0107       float nphot_long = 0;
0108       float nphot_short = 0;
0109       // cycle over photons in shower -------------------------------------------
0110       for (int iph = 0; iph < nphot; ++iph) {
0111         if (fiberId[iph] == 1) {
0112           nphot_long++;
0113         } else {
0114           nphot_short++;
0115           z[iph] = -z[iph];
0116         }
0117
0118         HFShowerPhoton::Point pos(x[iph], y[iph], z[iph]);
0119         HFShowerPhoton aPhoton(pos, t[iph], lambda[iph]);
0120         if (particle == "electron") {
0121           emColl.push_back(aPhoton);
0122         } else {
0123           hadColl.push_back(aPhoton);
0124         }
0125       }
0126       // end of cycle over photons in shower -------------------------------------------
0127
0128       if (particle == "electron") {
0129         LibTree->SetEntries(nem + 1);
0130         emBranch->Fill();
0131         nem++;
0132         emColl.clear();
0133       } else {
0134         LibTree->SetEntries(nhad + 1);
0135         nhad++;
0136         hadBranch->Fill();
0137         hadColl.clear();
0138       }
0139     }
0140     // end of cycle over showers ====================================================
0141     theFile->Close();
0142   }
0143   // end of cycle over files ++++++++++++++++++++++++++++++++++++++++++++++++++++
0144 }
0145
0146 void HcalForwardLibWriter::beginJob() {}
0147
0148 void HcalForwardLibWriter::endJob() {
0149   edm::Service<TFileService> fs;
0150   fs->file().cd();
0151   LibTree->Write();
0152   LibTree->Print();
0153 }
0154
0155 int HcalForwardLibWriter::readUserData(void) {
0156   std::ifstream input(theDataFile.c_str());
0157   if (input.fail()) {
0158     return 0;
0159   }
0160   std::string theFileName, thePID;
0161   int mom;
0162   int k = 0;
0163   while (!input.eof()) {
0164     input >> theFileName >> thePID >> mom;
0165     if (!input.fail()) {
0166       FileHandle aFile;
0167       aFile.name = theFileName;
0168       aFile.id = thePID;
0169       aFile.momentum = mom;
0170       theFileHandle.push_back(aFile);
0171       ++k;
0172     } else {
0173       input.clear();
0174     }
0175     input.ignore(999, '\n');
0176   }
0177   return k;
0178 }
0179
0180 //define this as a plug-in
0181 DEFINE_FWK_MODULE(HcalForwardLibWriter);