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     // cycle over showers ====================================================
0094     for (int iev = 0; iev < nentries; iev++) {
0095       if (primZ < 990.)
0096         continue;  // exclude showers with interactions in front of HF (1m of air)
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       // cycle over photons in shower -------------------------------------------
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       // end of cycle over photons in shower -------------------------------------------
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     // end of cycle over showers ====================================================
0139     theFile->Close();
0140   }
0141   // end of cycle over files ++++++++++++++++++++++++++++++++++++++++++++++++++++
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 //define this as a plug-in
0179 DEFINE_FWK_MODULE(HcalForwardLibWriter);