CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CaloMeanResponse.cc
Go to the documentation of this file.
3 
4 #include "CLHEP/Units/GlobalSystemOfUnits.h"
5 
6 #include <iostream>
7 #include <fstream>
8 #include <iomanip>
9 
10 //#define DebugLog
11 
13 
14  edm::ParameterSet m_p = p.getParameter<edm::ParameterSet>("CaloResponse");
15  edm::FileInPath fp = m_p.getParameter<edm::FileInPath>("ResponseFile");
16  std::string fName = fp.fullPath();
17  useTable = m_p.getParameter<bool>("UseResponseTable");
18  scale = m_p.getParameter<double>("ResponseScale");
19  edm::LogInfo("CaloSim") << "CaloMeanResponse initialized with scale "
20  << scale << " and use Table " << useTable
21  << " from file " << fName;
22 
23  readResponse (fName);
24 }
25 
27 
28 double CaloMeanResponse::getWeight(int genPID, double genP) {
29  double weight = 1;
30  bool found=false;
31  for (unsigned int i=0; i<pionTypes.size(); i++) {
32  if (genPID == pionTypes[i]) {
33  found = true;
34  break;
35  }
36  }
37  if (found) {
38  weight = scale;
39  if (useTable) {
40  if (piLast >= 0) weight *= pionTable[piLast];
41  for (unsigned int i=0; i<pionTable.size(); i++) {
42  if (genP < pionMomentum[i]) {
43  if (i == 0) weight = scale*pionTable[i];
44  else weight = scale*((pionTable[i-1]*(pionMomentum[i]-genP)+
45  pionTable[i]*(genP-pionMomentum[i-1]))/
46  (pionMomentum[i]-pionMomentum[i-1]));
47  break;
48  }
49  }
50  }
51 #ifdef DebugLog
52  LogDebug("CaloSim") << "CaloMeanResponse::getWeight PID " << genPID
53  << " uses pion list and gets weight " << weight
54  << " for momentum " << genP/GeV << " GeV/c";
55 #endif
56  } else {
57  for (unsigned int i=0; i<protonTypes.size(); i++) {
58  if (genPID == protonTypes[i]) {
59  found = true;
60  break;
61  }
62  }
63  if (found) {
64  weight = scale;
65  if (useTable) {
66  if (pLast >= 0) weight *= protonTable[pLast];
67  for (unsigned int i=0; i<protonTable.size(); i++) {
68  if (genP < protonMomentum[i]) {
69  if (i == 0) weight = scale*protonTable[i];
70  else weight = scale*((protonTable[i-1]*(protonMomentum[i]-genP)+
71  protonTable[i]*(genP-protonMomentum[i-1]))/
72  (protonMomentum[i]-protonMomentum[i-1]));
73  break;
74  }
75  }
76  }
77 #ifdef DebugLog
78  LogDebug("CaloSim") << "CaloMeanResponse::getWeight PID " << genPID
79  << " uses proton list and gets weight " << weight
80  << " for momentum " << genP/GeV << " GeV/c";
81 #endif
82  } else {
83 #ifdef DebugLog
84  LogDebug("CaloSim") << "CaloMeanResponse::getWeight PID " << genPID
85  << " is not in either lists and weight " << weight;
86 #endif
87  }
88  }
89  return weight;
90 }
91 
92 void CaloMeanResponse::readResponse (std::string fName) {
93 
94  std::ifstream infile;
95  infile.open(fName.c_str(), std::ios::in);
96 
97  if (infile) {
98  int nene, npart, pid;
99  double ene, responseData, responseMC, ratio;
100 
101  // First read the pion data
102  infile >> nene >> npart;
103  for (int i=0; i<npart; i++) {
104  infile >> pid;
105  pionTypes.push_back(pid);
106  }
107  for (int i=0; i<nene; i++) {
108  infile >> ene >> responseData >> responseMC;
109  if (responseMC > 0) ratio = responseData/responseMC;
110  else ratio = 1;
111  pionMomentum.push_back(ene*GeV);
112  pionTable.push_back(ratio);
113  }
114 
115  // Then read the proton data
116  infile >> nene >> npart;
117  for (int i=0; i<npart; i++) {
118  infile >> pid;
119  protonTypes.push_back(pid);
120  }
121  for (int i=0; i<nene; i++) {
122  infile >> ene >> responseData >> responseMC;
123  if (responseMC > 0) ratio = responseData/responseMC;
124  else ratio = 1;
125  protonMomentum.push_back(ene*GeV);
126  protonTable.push_back(ratio);
127  }
128  infile.close();
129  }
130 
131  piLast = (int)(pionTable.size()) - 1;
132  pLast = (int)(protonTable.size()) - 1;
133 #ifdef DebugLog
134  LogDebug("CaloSim") << "CaloMeanResponse::readResponse finds "
135  << pionTypes.size() << " particles to use pion response"
136  << " map with a table of " << pionTable.size()
137  << " data points " << piLast;
138  for (unsigned int i=0; i<pionTypes.size(); i++)
139  LogDebug("CaloSim") << "Particle ID[" << i << "] = " << pionTypes[i];
140  for (unsigned int i=0; i<pionTable.size(); i++)
141  LogDebug("CaloSim") << "Momentum[" << i << "] (" << pionMomentum[i]/GeV
142  << " GeV/c) --> " << pionTable[i];
143  LogDebug("CaloSim") << "CaloMeanResponse::readResponse finds "
144  << protonTypes.size() << " particles to use proton "
145  << "response map with a table of " << protonTable.size()
146  << " data points " << pLast;
147  for (unsigned int i=0; i<protonTypes.size(); i++)
148  LogDebug("CaloSim") << "Particle ID[" << i << "] = " << protonTypes[i];
149  for (unsigned int i=0; i<protonTable.size(); i++)
150  LogDebug("CaloSim") << "Momentum[" << i << "] (" << protonMomentum[i]/GeV
151  << " GeV/c) --> " << protonTable[i];
152 #endif
153 }
#define LogDebug(id)
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
void readResponse(std::string fName)
std::vector< int > protonTypes
double npart
Definition: HydjetWrapper.h:45
virtual ~CaloMeanResponse()
double getWeight(int genPID, double genP)
std::vector< double > protonMomentum
std::vector< double > pionTable
std::vector< int > pionTypes
CaloMeanResponse(edm::ParameterSet const &p)
list infile
Definition: EdgesToViz.py:90
std::vector< double > protonTable
std::vector< double > pionMomentum
std::string fullPath() const
Definition: FileInPath.cc:171