Back to home page

Project CMSSW displayed by LXR

 
 

    


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 //#define EDM_ML_DEBUG
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     // First read the pion data
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     // Then read the proton data
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 }