CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
PuppiAlgo.cc
Go to the documentation of this file.
3 #include "Math/QuantFuncMathCore.h"
4 #include "Math/SpecFuncMathCore.h"
5 #include "Math/ProbFunc.h"
6 #include "TMath.h"
7 
9  fEtaMin = iConfig.getParameter<std::vector<double>>("etaMin");
10  fEtaMax = iConfig.getParameter<std::vector<double>>("etaMax");
11  fPtMin = iConfig.getParameter<std::vector<double>>("ptMin");
12  fNeutralPtMin = iConfig.getParameter<std::vector<double>>("MinNeutralPt"); // Weighted Neutral Pt Cut
13  fNeutralPtSlope = iConfig.getParameter<std::vector<double>>("MinNeutralPtSlope"); // Slope vs #pv
14  fRMSEtaSF = iConfig.getParameter<std::vector<double>>("RMSEtaSF");
15  fMedEtaSF = iConfig.getParameter<std::vector<double>>("MedEtaSF");
16  fEtaMaxExtrap = iConfig.getParameter<double>("EtaMaxExtrap");
17 
18  std::vector<edm::ParameterSet> lAlgos = iConfig.getParameter<std::vector<edm::ParameterSet>>("puppiAlgos");
19  fNAlgos = lAlgos.size();
20  //Uber Configurable Puppi
21  std::vector<double> tmprms;
22  std::vector<double> tmpmed;
23 
24  for (unsigned int i0 = 0; i0 < lAlgos.size(); i0++) {
25  int pAlgoId = lAlgos[i0].getParameter<int>("algoId");
26  bool pCharged = lAlgos[i0].getParameter<bool>("useCharged");
27  bool pWeight0 = lAlgos[i0].getParameter<bool>("applyLowPUCorr");
28  int pComb = lAlgos[i0].getParameter<int>("combOpt"); // 0=> add in chi2/1=>Multiply p-values
29  double pConeSize = lAlgos[i0].getParameter<double>("cone"); // Min Pt when computing pt and rms
30  double pRMSPtMin = lAlgos[i0].getParameter<double>("rmsPtMin"); // Min Pt when computing pt and rms
31  double pRMSSF = lAlgos[i0].getParameter<double>("rmsScaleFactor"); // Additional Tuning parameter for Jokers
32  fAlgoId.push_back(pAlgoId);
33  fCharged.push_back(pCharged);
34  fAdjust.push_back(pWeight0);
35  fCombId.push_back(pComb);
36  fConeSize.push_back(pConeSize);
37  fRMSPtMin.push_back(pRMSPtMin);
38  fRMSScaleFactor.push_back(pRMSSF);
39  double pRMS = 0;
40  double pMed = 0;
41  double pMean = 0;
42  int pNCount = 0;
43  fRMS.push_back(pRMS);
44  fMedian.push_back(pMed);
45  fMean.push_back(pMean);
46  fNCount.push_back(pNCount);
47 
48  tmprms.clear();
49  tmpmed.clear();
50  for (unsigned int j0 = 0; j0 < fEtaMin.size(); j0++) {
51  tmprms.push_back(pRMS);
52  tmpmed.push_back(pMed);
53  }
54  fRMS_perEta.push_back(tmprms);
55  fMedian_perEta.push_back(tmpmed);
56  }
57 
58  cur_PtMin = -99.;
59  cur_NeutralPtMin = -99.;
60  cur_NeutralPtSlope = -99.;
61  cur_RMS = -99.;
62  cur_Med = -99.;
63 }
65  fPups.clear();
66  fPupsPV.clear();
67 }
69  fPups.clear();
70  fPupsPV.clear();
71  for (unsigned int i0 = 0; i0 < fNAlgos; i0++) {
72  fMedian[i0] = 0;
73  fRMS[i0] = 0;
74  fMean[i0] = 0;
75  fNCount[i0] = 0;
76  }
77 }
78 
79 void PuppiAlgo::fixAlgoEtaBin(int i_eta) {
80  cur_PtMin = fPtMin[i_eta];
83  cur_RMS = fRMS_perEta[0][i_eta]; // 0 is number of algos within this eta bin
84  cur_Med = fMedian_perEta[0][i_eta]; // 0 is number of algos within this eta bin
85 }
86 
87 void PuppiAlgo::add(const PuppiCandidate &iParticle, const double &iVal, const unsigned int iAlgo) {
88  if (iParticle.pt < fRMSPtMin[iAlgo])
89  return;
90  // Change from SRR : Previously used fastjet::PseudoJet::user_index to decide the particle type.
91  // In CMSSW we use the user_index to specify the index in the input collection, so I invented
92  // a new mechanism using the fastjet UserInfo functionality. Of course, it's still just an integer
93  // but that interface could be changed (or augmented) if desired / needed.
94  int puppi_id = iParticle.id;
95  if (puppi_id == std::numeric_limits<int>::lowest()) {
96  throw cms::Exception("PuppiRegisterNotSet") << "The puppi register is not set. This must be set before use.\n";
97  }
98 
99  // added by Nhan -- for all eta regions, compute mean/RMS from the central charged PU
100  if ((std::abs(iParticle.eta) < fEtaMaxExtrap) && (puppi_id == 2)) {
101  fPups.push_back(iVal);
102  fNCount[iAlgo]++;
103  }
104  // for the low PU case, correction. for checking that the PU-only median will be below the PV particles
105  if (std::abs(iParticle.eta) < fEtaMaxExtrap && (puppi_id == 1))
106  fPupsPV.push_back(iVal);
107 }
108 
111 //NHAN'S VERSION
112 void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo) {
113  //std::cout << "fNCount[iAlgo] = " << fNCount[iAlgo] << std::endl;
114  if (iAlgo >= fNAlgos)
115  return;
116  if (fNCount[iAlgo] == 0)
117  return;
118 
119  // sort alphas
120  int lNBefore = 0;
121  for (unsigned int i0 = 0; i0 < iAlgo; i0++)
122  lNBefore += fNCount[i0];
123  std::sort(fPups.begin() + lNBefore, fPups.begin() + lNBefore + fNCount[iAlgo]);
124 
125  // in case you have alphas == 0
126  int lNum0 = 0;
127  for (int i0 = lNBefore; i0 < lNBefore + fNCount[iAlgo]; i0++) {
128  if (fPups[i0] == 0)
129  lNum0 = i0 - lNBefore;
130  }
131 
132  // comput median, removed lCorr for now
133  int lNHalfway = lNBefore + lNum0 + int(double(fNCount[iAlgo] - lNum0) * 0.50);
134  fMedian[iAlgo] = fPups[lNHalfway];
135  double lMed = fMedian[iAlgo]; //Just to make the readability easier
136 
137  int lNRMS = 0;
138  for (int i0 = lNBefore; i0 < lNBefore + fNCount[iAlgo]; i0++) {
139  fMean[iAlgo] += fPups[i0];
140  if (fPups[i0] == 0)
141  continue;
142  // if(!fCharged[iAlgo] && fAdjust[iAlgo] && fPups[i0] > lMed) continue;
143  if (fAdjust[iAlgo] && fPups[i0] > lMed)
144  continue;
145  lNRMS++;
146  fRMS[iAlgo] += (fPups[i0] - lMed) * (fPups[i0] - lMed);
147  }
148  fMean[iAlgo] /= fNCount[iAlgo];
149  if (lNRMS > 0)
150  fRMS[iAlgo] /= lNRMS;
151  if (fRMS[iAlgo] == 0)
152  fRMS[iAlgo] = 1e-5;
153  // here is the raw RMS
154  fRMS[iAlgo] = sqrt(fRMS[iAlgo]);
155 
156  // some ways to do corrections to fRMS and fMedian
157  fRMS[iAlgo] *= fRMSScaleFactor[iAlgo];
158 
159  if (fAdjust[iAlgo]) {
160  //Adjust the p-value to correspond to the median
161  int lNPV = 0;
162  for (unsigned int i0 = 0; i0 < fPupsPV.size(); i0++)
163  if (fPupsPV[i0] <= lMed)
164  lNPV++;
165  double lAdjust = double(lNPV) / double(lNPV + 0.5 * fNCount[iAlgo]);
166  if (lAdjust > 0) {
167  fMedian[iAlgo] -= sqrt(ROOT::Math::chisquared_quantile(lAdjust, 1.) * fRMS[iAlgo]);
168  fRMS[iAlgo] -= sqrt(ROOT::Math::chisquared_quantile(lAdjust, 1.) * fRMS[iAlgo]);
169  }
170  }
171 
172  // fRMS_perEta[iAlgo] *= cur_RMSEtaSF;
173  // fMedian_perEta[iAlgo] *= cur_MedEtaSF;
174 
175  for (unsigned int j0 = 0; j0 < fEtaMin.size(); j0++) {
176  fRMS_perEta[iAlgo][j0] = fRMS[iAlgo] * fRMSEtaSF[j0];
177  fMedian_perEta[iAlgo][j0] = fMedian[iAlgo] * fMedEtaSF[j0];
178  }
179 }
181 
182 //This code is probably a bit confusing
183 double PuppiAlgo::compute(std::vector<double> const &iVals, double iChi2) const {
184  if (fAlgoId[0] == -1)
185  return 1;
186  double lVal = 0.;
187  double lPVal = 1.;
188  int lNDOF = 0;
189  for (unsigned int i0 = 0; i0 < fNAlgos; i0++) {
190  if (fNCount[i0] == 0)
191  return 1.; //in the NoPU case return 1.
192  if (fCombId[i0] == 1 && i0 > 0) { //Compute the previous p-value so that p-values can be multiplieed
193  double pPVal = ROOT::Math::chisquared_cdf(lVal, lNDOF);
194  lPVal *= pPVal;
195  lNDOF = 0;
196  lVal = 0;
197  }
198  double pVal = iVals[i0];
199  //Special Check for any algo with log(0)
200  if (fAlgoId[i0] == 0 && iVals[i0] == 0)
201  pVal = cur_Med;
202  if (fAlgoId[i0] == 3 && iVals[i0] == 0)
203  pVal = cur_Med;
204  if (fAlgoId[i0] == 5 && iVals[i0] == 0)
205  pVal = cur_Med;
206  lVal += (pVal - cur_Med) * (fabs(pVal - cur_Med)) / cur_RMS / cur_RMS;
207  lNDOF++;
208  if (i0 == 0 && iChi2 != 0)
209  lNDOF++; //Add external Chi2 to first element
210  if (i0 == 0 && iChi2 != 0)
211  lVal += iChi2; //Add external Chi2 to first element
212  }
213  //Top it off with the last calc
214  lPVal *= ROOT::Math::chisquared_cdf(lVal, lNDOF);
215  return lPVal;
216 }
217 // ------------------------------------------------------------------------------------------
219  edm::ParameterSetDescription puppialgos;
220  puppialgos.add<int>("algoId", 5);
221  puppialgos.add<bool>("useCharged", false);
222  puppialgos.add<bool>("applyLowPUCorr", false);
223  puppialgos.add<int>("combOpt", 5);
224  puppialgos.add<double>("cone", .4);
225  puppialgos.add<double>("rmsPtMin", .1);
226  puppialgos.add<double>("rmsScaleFactor", 1.0);
227  std::vector<edm::ParameterSet> VPSetPuppiAlgos;
228  edm::ParameterSet puppiset;
229  puppiset.addParameter<int>("algoId", 5);
230  puppiset.addParameter<bool>("useCharged", false);
231  puppiset.addParameter<bool>("applyLowPUCorr", false);
232  puppiset.addParameter<int>("combOpt", 5);
233  puppiset.addParameter<double>("cone", .4);
234  puppiset.addParameter<double>("rmsPtMin", .1);
235  puppiset.addParameter<double>("rmsScaleFactor", 1.0);
236  VPSetPuppiAlgos.push_back(puppiset);
237 
239  algos.addVPSet("puppiAlgos", puppialgos, VPSetPuppiAlgos);
240  std::vector<edm::ParameterSet> VPSetAlgos;
241  edm::ParameterSet algosset;
242  algos.add<std::vector<double>>("etaMin", {0.});
243  algos.add<std::vector<double>>("etaMax", {2.5});
244  algos.add<std::vector<double>>("ptMin", {0.});
245  algos.add<std::vector<double>>("MinNeutralPt", {0.2});
246  algos.add<std::vector<double>>("MinNeutralPtSlope", {0.015});
247  algos.add<std::vector<double>>("RMSEtaSF", {1.0});
248  algos.add<std::vector<double>>("MedEtaSF", {1.0});
249  algos.add<double>("EtaMaxExtrap", 2.0);
250  algosset.addParameter<std::vector<double>>("etaMin", {0.});
251  algosset.addParameter<std::vector<double>>("etaMax", {2.5});
252  algosset.addParameter<std::vector<double>>("ptMin", {0.});
253  algosset.addParameter<std::vector<double>>("MinNeutralPt", {0.2});
254  algosset.addParameter<std::vector<double>>("MinNeutralPtSlope", {0.015});
255  algosset.addParameter<std::vector<double>>("RMSEtaSF", {1.0});
256  algosset.addParameter<std::vector<double>>("MedEtaSF", {1.0});
257  algosset.addParameter<double>("EtaMaxExtrap", 2.0);
258  algosset.addParameter<std::vector<edm::ParameterSet>>("puppiAlgos", VPSetPuppiAlgos);
259  VPSetAlgos.push_back(algosset);
260  desc.addVPSet("algos", algos, VPSetAlgos);
261 }
void add(const PuppiCandidate &iParticle, const double &iVal, const unsigned int iAlgo)
Definition: PuppiAlgo.cc:87
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
void computeMedRMS(const unsigned int &iAlgo)
Definition: PuppiAlgo.cc:112
double cur_PtMin
Definition: PuppiAlgo.h:52
std::vector< double > fEtaMin
Definition: PuppiAlgo.h:43
static void fillDescriptionsPuppiAlgo(edm::ParameterSetDescription &desc)
Definition: PuppiAlgo.cc:218
std::vector< bool > fAdjust
Definition: PuppiAlgo.h:67
PuppiAlgo(edm::ParameterSet &iConfig)
Definition: PuppiAlgo.cc:8
double cur_NeutralPtSlope
Definition: PuppiAlgo.h:54
std::vector< double > fConeSize
Definition: PuppiAlgo.h:69
void reset()
Definition: PuppiAlgo.cc:68
double cur_Med
Definition: PuppiAlgo.h:56
std::vector< std::vector< double > > fMedian_perEta
Definition: PuppiAlgo.h:61
double compute(std::vector< double > const &iVals, double iChi2) const
Definition: PuppiAlgo.cc:183
T sqrt(T t)
Definition: SSEVec.h:19
~PuppiAlgo()
Definition: PuppiAlgo.cc:64
std::vector< double > fRMS
Definition: PuppiAlgo.h:58
void addParameter(std::string const &name, T const &value)
Definition: ParameterSet.h:135
std::vector< double > fNeutralPtMin
Definition: PuppiAlgo.h:45
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< float > fPupsPV
Definition: PuppiAlgo.h:64
ParameterDescriptionBase * add(U const &iLabel, T const &value)
double cur_NeutralPtMin
Definition: PuppiAlgo.h:53
std::vector< float > fPups
Definition: PuppiAlgo.h:63
std::vector< double > fMean
Definition: PuppiAlgo.h:72
std::vector< int > fAlgoId
Definition: PuppiAlgo.h:65
std::vector< double > fRMSPtMin
Definition: PuppiAlgo.h:70
std::vector< int > fCombId
Definition: PuppiAlgo.h:68
std::vector< double > fRMSScaleFactor
Definition: PuppiAlgo.h:71
std::vector< double > fPtMin
Definition: PuppiAlgo.h:44
std::vector< double > fRMSEtaSF
Definition: PuppiAlgo.h:48
std::vector< double > fNeutralPtSlope
Definition: PuppiAlgo.h:46
std::vector< double > fMedEtaSF
Definition: PuppiAlgo.h:49
std::vector< double > fEtaMax
Definition: PuppiAlgo.h:42
double fEtaMaxExtrap
Definition: PuppiAlgo.h:50
double cur_RMS
Definition: PuppiAlgo.h:55
void fixAlgoEtaBin(int i_eta)
Definition: PuppiAlgo.cc:79
std::vector< bool > fCharged
Definition: PuppiAlgo.h:66
std::vector< std::vector< double > > fRMS_perEta
Definition: PuppiAlgo.h:60
unsigned int fNAlgos
Definition: PuppiAlgo.h:41
std::vector< int > fNCount
Definition: PuppiAlgo.h:73
std::vector< double > fMedian
Definition: PuppiAlgo.h:59