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