Fork me on GitHub

source: git/external/PUPPI/PuppiAlgo.cc@ 0882413

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 0882413 was 659c7b6, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago

fix PUPPI includes

  • Property mode set to 100644
File size: 6.8 KB
Line 
1#include "PuppiAlgo.hh"
2#include "Math/QuantFuncMathCore.h"
3#include "Math/SpecFuncMathCore.h"
4#include "Math/ProbFunc.h"
5#include "TMath.h"
6
7PuppiAlgo::PuppiAlgo(AlgoObj &iAlgo) {
8 fEtaMin = iAlgo.etaMin;
9 fEtaMax = iAlgo.etaMax;
10 fPtMin = iAlgo.ptMin;
11 fNeutralPtMin = iAlgo.minNeutralPt;
12 fNeutralPtSlope = iAlgo.minNeutralPtSlope;
13 fRMSEtaSF = iAlgo.rmsEtaSF;
14 fMedEtaSF = iAlgo.medEtaSF;
15 fEtaMaxExtrap = iAlgo.etaMaxExtrap;
16 std::vector<AlgoSubObj> lAlgos = iAlgo.subAlgos;
17 fNAlgos = lAlgos.size();
18 //std::cout << "==> " << fEtaMin << " - " << fEtaMax << " - " << fPtMin << " - " << fNeutralPtMin << " - " << fNeutralPtSlope << " - " << fRMSEtaSF << " - " << std::endl;
19 for(unsigned int i0 = 0; i0 < lAlgos.size(); i0++) {
20 int pAlgoId = lAlgos[i0].metricId;
21 bool pCharged = lAlgos[i0].useCharged;
22 bool pWeight0 = lAlgos[i0].applyLowPUCorr;
23 int pComb = lAlgos[i0].combId;
24 double pConeSize = lAlgos[i0].coneSize;
25 double pRMSPtMin = lAlgos[i0].rmsPtMin;
26 double pRMSSF = lAlgos[i0].rmsScaleFactor;
27 //std::cout << "Algo==> " << i0 << " - " << pAlgoId << " - " << pCharged << " - " << pWeight0 << " - " << pComb << " - " << pConeSize << " - " << pRMSPtMin << " - " << std::endl;
28 fAlgoId .push_back(pAlgoId);
29 fCharged .push_back(pCharged);
30 fAdjust .push_back(pWeight0);
31 fCombId .push_back(pComb);
32 fConeSize .push_back(pConeSize);
33 fRMSPtMin .push_back(pRMSPtMin);
34 fRMSScaleFactor.push_back(pRMSSF);
35 double pRMS = 0;
36 double pMed = 0;
37 double pMean = 0;
38 int pNCount = 0;
39 fRMS .push_back(pRMS);
40 fMedian.push_back(pMed);
41 fMean .push_back(pMean);
42 fNCount.push_back(pNCount);
43 /*
44 tmprms .clear();
45 tmpmed .clear();
46 for (unsigned int j0 = 0; j0 < fEtaMin.size(); j0++){
47 tmprms.push_back(pRMS);
48 tmpmed.push_back(pMed);
49 }
50 fRMS_perEta.push_back(tmprms);
51 fMedian_perEta.push_back(tmpmed);
52 */
53 }
54}
55PuppiAlgo::~PuppiAlgo() {
56 fPups .clear();
57 fPupsPV.clear();
58}
59void PuppiAlgo::reset() {
60 fPups .clear();
61 fPupsPV.clear();
62 for(unsigned int i0 = 0; i0 < fNAlgos; i0++) {
63 fMedian[i0] = 0;
64 fRMS [i0] = 0;
65 fMean [i0] = 0;
66 fNCount[i0] = 0;
67 }
68}
69void PuppiAlgo::add(const fastjet::PseudoJet &iParticle,const double &iVal,const unsigned int iAlgo) {
70 if(iParticle.pt() < fRMSPtMin[iAlgo]) return;
71 if(fCharged[iAlgo] && fabs(iParticle.user_index()) < 1) return;
72 if(fCharged[iAlgo] && (fabs(iParticle.user_index()) >=1 && fabs(iParticle.user_index()) <=2)) fPupsPV.push_back(iVal);
73 if(fCharged[iAlgo] && fabs(iParticle.user_index()) < 3) return;
74 fPups.push_back(iVal);
75 fNCount[iAlgo]++;
76 /* eta extrap (cmssw version)
77 if ((std::abs(iParticle.eta()) < fEtaMaxExtrap) && (std::abs(puppi_register) >= 3)){
78 fPups.push_back(iVal);
79 // fPupsPV.push_back(iVal);
80 fNCount[iAlgo]++;
81 }
82 // for the low PU case, correction. for checking that the PU-only median will be below the PV particles
83 if(std::abs(iParticle.eta()) < fEtaMaxExtrap && (std::abs(puppi_register) >=1 && std::abs(puppi_register) <=2)) fPupsPV.push_back(iVal);
84 */
85}
86void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo,const double &iPVFrac) {
87 if(iAlgo >= fNAlgos ) return;
88 if(fNCount[iAlgo] == 0) return;
89 int lNBefore = 0;
90 for(unsigned int i0 = 0; i0 < iAlgo; i0++) lNBefore += fNCount[i0];
91 std::sort(fPups.begin()+lNBefore,fPups.begin()+lNBefore+fNCount[iAlgo]);
92 int lNum0 = 0;
93 for(int i0 = lNBefore; i0 < lNBefore+fNCount[iAlgo]; i0++) {
94 if(fPups[i0] == 0) lNum0 = i0-lNBefore;
95 }
96 //lNum0 = 0;
97 int lNHalfway = lNBefore + lNum0 + int( double( fNCount[iAlgo]-lNum0 )*0.50);
98 fMedian[iAlgo] = fPups[lNHalfway];
99 double lMed = fMedian[iAlgo]; //Just to make the readability easier
100
101 int lNRMS = 0;
102 for(int i0 = lNBefore; i0 < lNBefore+fNCount[iAlgo]; i0++) {
103 fMean[iAlgo] += fPups[i0];
104 if(fPups[i0] == 0) continue;
105 if(fAdjust[iAlgo] && fPups[i0] > lMed) continue;
106 lNRMS++;
107 fRMS [iAlgo] += (fPups[i0]-lMed)*(fPups[i0]-lMed);
108 }
109 fMean[iAlgo]/=fNCount[iAlgo];
110 if(lNRMS > 0) fRMS [iAlgo]/=lNRMS;
111 if(fRMS[iAlgo] == 0) fRMS[iAlgo] = 1e-5;
112
113 fRMS [iAlgo] = sqrt(fRMS[iAlgo]);
114 fRMS [iAlgo] *= fRMSScaleFactor[iAlgo];
115 if(fAdjust[iAlgo]){
116 //Adjust the p-value to correspond to the median
117 std::sort(fPupsPV.begin(),fPupsPV.end());
118 int lNPV = 0;
119 for(unsigned int i0 = 0; i0 < fPupsPV.size(); i0++) if(fPupsPV[i0] <= lMed ) lNPV++;
120 double lAdjust = double(lNPV)/double(lNPV+0.5*fNCount[iAlgo]);
121 if(lAdjust > 0) {
122 fMedian[iAlgo] -= sqrt(ROOT::Math::chisquared_quantile(lAdjust,1.)*fRMS[iAlgo]);
123 fRMS[iAlgo] -= sqrt(ROOT::Math::chisquared_quantile(lAdjust,1.)*fRMS[iAlgo]);
124 }
125 }
126 /* eta extrapolated version
127 for (unsigned int j0 = 0; j0 < fEtaMin.size(); j0++){
128 fRMS_perEta[iAlgo][j0] = fRMS[iAlgo]*fRMSEtaSF[j0];
129 fMedian_perEta[iAlgo][j0] = fMedian[iAlgo]*fMedEtaSF[j0];
130 }
131 */
132}
133//This code is probably a bit confusing
134double PuppiAlgo::compute(std::vector<double> &iVals,double iChi2) {
135 if(fAlgoId[0] == -1) return 1;
136 double lVal = 0.;
137 double lPVal = 1.;
138 int lNDOF = 0;
139 for(unsigned int i0 = 0; i0 < fNAlgos; i0++) {
140 if(fNCount[i0] == 0) return 1.; //in the NoPU case return 1.
141 if(fCombId[i0] == 1 && i0 > 0) { //Compute the previous p-value so that p-values can be multiplieed
142 double pPVal = ROOT::Math::chisquared_cdf(lVal,lNDOF);
143 lPVal *= pPVal;
144 lNDOF = 0;
145 lVal = 0;
146 }
147 double pVal = iVals[i0];
148 //Special Check for any algo with log(0)
149 if(fAlgoId[i0] == 0 && iVals[i0] == 0) pVal = fMedian[i0];
150 if(fAlgoId[i0] == 3 && iVals[i0] == 0) pVal = fMedian[i0];
151 if(fAlgoId[i0] == 5 && iVals[i0] == 0) pVal = fMedian[i0];
152 lVal += (pVal-fMedian[i0])*(fabs(pVal-fMedian[i0]))/fRMS[i0]/fRMS[i0];
153 lNDOF++;
154 if(i0 == 0 && iChi2 != 0) lNDOF++; //Add external Chi2 to first element
155 if(i0 == 0 && iChi2 != 0) lVal+=iChi2; //Add external Chi2 to first element
156 }
157 //Top it off with the last calc
158 lPVal *= ROOT::Math::chisquared_cdf(lVal,lNDOF);
159 return lPVal;
160}
161double PuppiAlgo::neutralPt(int iNPV) {
162 return fNeutralPtMin + iNPV * fNeutralPtSlope;
163}
164int PuppiAlgo::numAlgos() {
165 return fNAlgos;
166}
167double PuppiAlgo::ptMin() {
168 return fPtMin;
169}
170double PuppiAlgo::etaMin() {
171 return fEtaMin;
172}
173double PuppiAlgo::etaMax() {
174 return fEtaMax;
175}
176int PuppiAlgo::algoId(const unsigned int &iAlgo) {
177 assert(iAlgo < fNAlgos);
178 return fAlgoId[iAlgo];
179}
180bool PuppiAlgo::isCharged(const unsigned int &iAlgo) {
181 assert(iAlgo < fNAlgos);
182 return fCharged[iAlgo];
183}
184double PuppiAlgo::coneSize(const unsigned int &iAlgo) {
185 assert(iAlgo < fNAlgos);
186 return fConeSize[iAlgo];
187}
Note: See TracBrowser for help on using the repository browser.