File indexing completed on 2024-12-20 03:14:12
0001 #include "SimG4CMS/Calo/interface/CaloMeanResponse.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003
0004 #include <CLHEP/Units/SystemOfUnits.h>
0005 using CLHEP::GeV;
0006
0007 #include <iostream>
0008 #include <fstream>
0009 #include <iomanip>
0010
0011
0012
0013 CaloMeanResponse::CaloMeanResponse(edm::ParameterSet const& p) {
0014 edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("CaloResponse");
0015 edm::FileInPath fp = m_p.getParameter<edm::FileInPath>("ResponseFile");
0016 const std::string& fName = fp.fullPath();
0017 useTable = m_p.getParameter<bool>("UseResponseTable");
0018 scale = m_p.getParameter<double>("ResponseScale");
0019 edm::LogVerbatim("CaloSim") << "CaloMeanResponse initialized with scale " << scale << " and use Table " << useTable
0020 << " from file " << fName;
0021 readResponse(fName);
0022 }
0023
0024 CaloMeanResponse::~CaloMeanResponse() {}
0025
0026 double CaloMeanResponse::getWeight(int genPID, double genP) {
0027 double weight = 1;
0028 bool found = false;
0029 for (unsigned int i = 0; i < pionTypes.size(); i++) {
0030 if (genPID == pionTypes[i]) {
0031 found = true;
0032 break;
0033 }
0034 }
0035 if (found) {
0036 weight = scale;
0037 if (useTable) {
0038 if (piLast >= 0)
0039 weight *= pionTable[piLast];
0040 for (unsigned int i = 0; i < pionTable.size(); i++) {
0041 if (genP < pionMomentum[i]) {
0042 if (i == 0)
0043 weight = scale * pionTable[i];
0044 else
0045 weight =
0046 scale * ((pionTable[i - 1] * (pionMomentum[i] - genP) + pionTable[i] * (genP - pionMomentum[i - 1])) /
0047 (pionMomentum[i] - pionMomentum[i - 1]));
0048 break;
0049 }
0050 }
0051 }
0052 #ifdef EDM_ML_DEBUG
0053 edm::LogVerbatim("CaloSim") << "CaloMeanResponse::getWeight PID " << genPID << " uses pion list and gets weight "
0054 << weight << " for momentum " << genP / GeV << " GeV/c";
0055 #endif
0056 } else {
0057 for (unsigned int i = 0; i < protonTypes.size(); i++) {
0058 if (genPID == protonTypes[i]) {
0059 found = true;
0060 break;
0061 }
0062 }
0063 if (found) {
0064 weight = scale;
0065 if (useTable) {
0066 if (pLast >= 0)
0067 weight *= protonTable[pLast];
0068 for (unsigned int i = 0; i < protonTable.size(); i++) {
0069 if (genP < protonMomentum[i]) {
0070 if (i == 0)
0071 weight = scale * protonTable[i];
0072 else
0073 weight =
0074 scale *
0075 ((protonTable[i - 1] * (protonMomentum[i] - genP) + protonTable[i] * (genP - protonMomentum[i - 1])) /
0076 (protonMomentum[i] - protonMomentum[i - 1]));
0077 break;
0078 }
0079 }
0080 }
0081 #ifdef EDM_ML_DEBUG
0082 edm::LogVerbatim("CaloSim") << "CaloMeanResponse::getWeight PID " << genPID << " uses proton list and gets "
0083 << "weight " << weight << " for momentum " << genP / GeV << " GeV/c";
0084 #endif
0085 } else {
0086 #ifdef EDM_ML_DEBUG
0087 edm::LogVerbatim("CaloSim") << "CaloMeanResponse::getWeight PID " << genPID << " is not in either lists and "
0088 << "weight " << weight;
0089 #endif
0090 }
0091 }
0092 return weight;
0093 }
0094
0095 void CaloMeanResponse::readResponse(std::string fName) {
0096 std::ifstream infile;
0097 infile.open(fName.c_str(), std::ios::in);
0098
0099 if (infile) {
0100 int nene, npart, pid;
0101 double ene, responseData, responseMC, ratio;
0102
0103
0104 infile >> nene >> npart;
0105 for (int i = 0; i < npart; i++) {
0106 infile >> pid;
0107 pionTypes.push_back(pid);
0108 }
0109 for (int i = 0; i < nene; i++) {
0110 infile >> ene >> responseData >> responseMC;
0111 if (responseMC > 0)
0112 ratio = responseData / responseMC;
0113 else
0114 ratio = 1;
0115 pionMomentum.push_back(ene * GeV);
0116 pionTable.push_back(ratio);
0117 }
0118
0119
0120 infile >> nene >> npart;
0121 for (int i = 0; i < npart; i++) {
0122 infile >> pid;
0123 protonTypes.push_back(pid);
0124 }
0125 for (int i = 0; i < nene; i++) {
0126 infile >> ene >> responseData >> responseMC;
0127 if (responseMC > 0)
0128 ratio = responseData / responseMC;
0129 else
0130 ratio = 1;
0131 protonMomentum.push_back(ene * GeV);
0132 protonTable.push_back(ratio);
0133 }
0134 infile.close();
0135 }
0136
0137 piLast = (int)(pionTable.size()) - 1;
0138 pLast = (int)(protonTable.size()) - 1;
0139 #ifdef EDM_ML_DEBUG
0140 edm::LogVerbatim("CaloSim") << "CaloMeanResponse::readResponse finds " << pionTypes.size()
0141 << " particles to use pion "
0142 << "response map with a table of " << pionTable.size() << " data points " << piLast;
0143 for (unsigned int i = 0; i < pionTypes.size(); i++)
0144 edm::LogVerbatim("CaloSim") << "Particle ID[" << i << "] = " << pionTypes[i];
0145 for (unsigned int i = 0; i < pionTable.size(); i++)
0146 edm::LogVerbatim("CaloSim") << "Momentum[" << i << "] (" << pionMomentum[i] / GeV << " GeV/c) --> " << pionTable[i];
0147 edm::LogVerbatim("CaloSim") << "CaloMeanResponse::readResponse finds " << protonTypes.size() << " particles to use "
0148 << "proton response map with a table of " << protonTable.size() << " data points "
0149 << pLast;
0150 for (unsigned int i = 0; i < protonTypes.size(); i++)
0151 edm::LogVerbatim("CaloSim") << "Particle ID[" << i << "] = " << protonTypes[i];
0152 for (unsigned int i = 0; i < protonTable.size(); i++)
0153 edm::LogVerbatim("CaloSim") << "Momentum[" << i << "] (" << protonMomentum[i] / GeV << " GeV/c) --> "
0154 << protonTable[i];
0155 #endif
0156 }
