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