| 1 | #include "PuppiAlgo.hh"
|
|---|
| 2 | #include "Math/QuantFuncMathCore.h"
|
|---|
| 3 | #include "Math/SpecFuncMathCore.h"
|
|---|
| 4 | #include "Math/ProbFunc.h"
|
|---|
| 5 | #include "TMath.h"
|
|---|
| 6 |
|
|---|
| 7 | PuppiAlgo::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 | }
|
|---|
| 55 | PuppiAlgo::~PuppiAlgo() {
|
|---|
| 56 | fPups .clear();
|
|---|
| 57 | fPupsPV.clear();
|
|---|
| 58 | }
|
|---|
| 59 | void 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 | }
|
|---|
| 69 | void 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 | }
|
|---|
| 86 | void 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
|
|---|
| 134 | double 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 | }
|
|---|
| 161 | double PuppiAlgo::neutralPt(int iNPV) {
|
|---|
| 162 | return fNeutralPtMin + iNPV * fNeutralPtSlope;
|
|---|
| 163 | }
|
|---|
| 164 | int PuppiAlgo::numAlgos() {
|
|---|
| 165 | return fNAlgos;
|
|---|
| 166 | }
|
|---|
| 167 | double PuppiAlgo::ptMin() {
|
|---|
| 168 | return fPtMin;
|
|---|
| 169 | }
|
|---|
| 170 | double PuppiAlgo::etaMin() {
|
|---|
| 171 | return fEtaMin;
|
|---|
| 172 | }
|
|---|
| 173 | double PuppiAlgo::etaMax() {
|
|---|
| 174 | return fEtaMax;
|
|---|
| 175 | }
|
|---|
| 176 | int PuppiAlgo::algoId(const unsigned int &iAlgo) {
|
|---|
| 177 | assert(iAlgo < fNAlgos);
|
|---|
| 178 | return fAlgoId[iAlgo];
|
|---|
| 179 | }
|
|---|
| 180 | bool PuppiAlgo::isCharged(const unsigned int &iAlgo) {
|
|---|
| 181 | assert(iAlgo < fNAlgos);
|
|---|
| 182 | return fCharged[iAlgo];
|
|---|
| 183 | }
|
|---|
| 184 | double PuppiAlgo::coneSize(const unsigned int &iAlgo) {
|
|---|
| 185 | assert(iAlgo < fNAlgos);
|
|---|
| 186 | return fConeSize[iAlgo];
|
|---|
| 187 | }
|
|---|
