Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:01:08 

0001 #include "CommonTools/PileupAlgos/interface/PuppiAlgo.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include "Math/QuantFuncMathCore.h"
0004 #include "Math/SpecFuncMathCore.h"
0005 #include "Math/ProbFunc.h"
0006 #include "TMath.h"
0007 
0008 PuppiAlgo::PuppiAlgo(edm::ParameterSet &iConfig) {
0009   fEtaMin = iConfig.getParameter<std::vector<double>>("etaMin");
0010   fEtaMax = iConfig.getParameter<std::vector<double>>("etaMax");
0011   fPtMin = iConfig.getParameter<std::vector<double>>("ptMin");
0012   fNeutralPtMin = iConfig.getParameter<std::vector<double>>("MinNeutralPt");         // Weighted Neutral Pt Cut
0013   fNeutralPtSlope = iConfig.getParameter<std::vector<double>>("MinNeutralPtSlope");  // Slope vs #pv
0014   fRMSEtaSF = iConfig.getParameter<std::vector<double>>("RMSEtaSF");
0015   fMedEtaSF = iConfig.getParameter<std::vector<double>>("MedEtaSF");
0016   fEtaMaxExtrap = iConfig.getParameter<double>("EtaMaxExtrap");
0017 
0018   std::vector<edm::ParameterSet> lAlgos = iConfig.getParameter<std::vector<edm::ParameterSet>>("puppiAlgos");
0019   fNAlgos = lAlgos.size();
0020   //Uber Configurable Puppi
0021   std::vector<double> tmprms;
0022   std::vector<double> tmpmed;
0023 
0024   for (unsigned int i0 = 0; i0 < lAlgos.size(); i0++) {
0025     int pAlgoId = lAlgos[i0].getParameter<int>("algoId");
0026     bool pCharged = lAlgos[i0].getParameter<bool>("useCharged");
0027     bool pWeight0 = lAlgos[i0].getParameter<bool>("applyLowPUCorr");
0028     int pComb = lAlgos[i0].getParameter<int>("combOpt");                // 0=> add in chi2/1=>Multiply p-values
0029     double pConeSize = lAlgos[i0].getParameter<double>("cone");         // Min Pt when computing pt and rms
0030     double pRMSPtMin = lAlgos[i0].getParameter<double>("rmsPtMin");     // Min Pt when computing pt and rms
0031     double pRMSSF = lAlgos[i0].getParameter<double>("rmsScaleFactor");  // Additional Tuning parameter for Jokers
0032     fAlgoId.push_back(pAlgoId);
0033     fCharged.push_back(pCharged);
0034     fAdjust.push_back(pWeight0);
0035     fCombId.push_back(pComb);
0036     fConeSize.push_back(pConeSize);
0037     fRMSPtMin.push_back(pRMSPtMin);
0038     fRMSScaleFactor.push_back(pRMSSF);
0039     double pRMS = 0;
0040     double pMed = 0;
0041     double pMean = 0;
0042     int pNCount = 0;
0043     fRMS.push_back(pRMS);
0044     fMedian.push_back(pMed);
0045     fMean.push_back(pMean);
0046     fNCount.push_back(pNCount);
0047 
0048     tmprms.clear();
0049     tmpmed.clear();
0050     for (unsigned int j0 = 0; j0 < fEtaMin.size(); j0++) {
0051       tmprms.push_back(pRMS);
0052       tmpmed.push_back(pMed);
0053     }
0054     fRMS_perEta.push_back(tmprms);
0055     fMedian_perEta.push_back(tmpmed);
0056   }
0057 
0058   cur_PtMin = -99.;
0059   cur_NeutralPtMin = -99.;
0060   cur_NeutralPtSlope = -99.;
0061   cur_RMS = -99.;
0062   cur_Med = -99.;
0063 }
0064 PuppiAlgo::~PuppiAlgo() {
0065   fPups.clear();
0066   fPupsPV.clear();
0067 }
0068 void PuppiAlgo::reset() {
0069   fPups.clear();
0070   fPupsPV.clear();
0071   for (unsigned int i0 = 0; i0 < fNAlgos; i0++) {
0072     fMedian[i0] = 0;
0073     fRMS[i0] = 0;
0074     fMean[i0] = 0;
0075     fNCount[i0] = 0;
0076   }
0077 }
0078 
0079 void PuppiAlgo::fixAlgoEtaBin(int i_eta) {
0080   cur_PtMin = fPtMin[i_eta];
0081   cur_NeutralPtMin = fNeutralPtMin[i_eta];
0082   cur_NeutralPtSlope = fNeutralPtSlope[i_eta];
0083   cur_RMS = fRMS_perEta[0][i_eta];     // 0 is number of algos within this eta bin
0084   cur_Med = fMedian_perEta[0][i_eta];  // 0 is number of algos within this eta bin
0085 }
0086 
0087 void PuppiAlgo::add(const PuppiCandidate &iParticle, const double &iVal, const unsigned int iAlgo) {
0088   if (iParticle.pt < fRMSPtMin[iAlgo])
0089     return;
0090   // Change from SRR : Previously used fastjet::PseudoJet::user_index to decide the particle type.
0091   // In CMSSW we use the user_index to specify the index in the input collection, so I invented
0092   // a new mechanism using the fastjet UserInfo functionality. Of course, it's still just an integer
0093   // but that interface could be changed (or augmented) if desired / needed.
0094   int puppi_id = iParticle.id;
0095   if (puppi_id == std::numeric_limits<int>::lowest()) {
0096     throw cms::Exception("PuppiRegisterNotSet") << "The puppi register is not set. This must be set before use.\n";
0097   }
0098 
0099   // added by Nhan -- for all eta regions, compute mean/RMS from the central charged PU
0100   if ((std::abs(iParticle.eta) < fEtaMaxExtrap) && (puppi_id == 2)) {
0101     fPups.push_back(iVal);
0102     fNCount[iAlgo]++;
0103   }
0104   // for the low PU case, correction.  for checking that the PU-only median will be below the PV particles
0105   if (std::abs(iParticle.eta) < fEtaMaxExtrap && (puppi_id == 1))
0106     fPupsPV.push_back(iVal);
0107 }
0108 
0109 ////////////////////////////////////////////////////////////////////////////////
0110 ////////////////////////////////////////////////////////////////////////////////
0111 //NHAN'S VERSION
0112 void PuppiAlgo::computeMedRMS(const unsigned int &iAlgo) {
0113   //std::cout << "fNCount[iAlgo] = " << fNCount[iAlgo] << std::endl;
0114   if (iAlgo >= fNAlgos)
0115     return;
0116   if (fNCount[iAlgo] == 0)
0117     return;
0118 
0119   // sort alphas
0120   int lNBefore = 0;
0121   for (unsigned int i0 = 0; i0 < iAlgo; i0++)
0122     lNBefore += fNCount[i0];
0123   std::sort(fPups.begin() + lNBefore, fPups.begin() + lNBefore + fNCount[iAlgo]);
0124 
0125   // in case you have alphas == 0
0126   int lNum0 = 0;
0127   for (int i0 = lNBefore; i0 < lNBefore + fNCount[iAlgo]; i0++) {
0128     if (fPups[i0] == 0)
0129       lNum0 = i0 - lNBefore;
0130   }
0131 
0132   // comput median, removed lCorr for now
0133   int lNHalfway = lNBefore + lNum0 + int(double(fNCount[iAlgo] - lNum0) * 0.50);
0134   fMedian[iAlgo] = fPups[lNHalfway];
0135   double lMed = fMedian[iAlgo];  //Just to make the readability easier
0136 
0137   int lNRMS = 0;
0138   for (int i0 = lNBefore; i0 < lNBefore + fNCount[iAlgo]; i0++) {
0139     fMean[iAlgo] += fPups[i0];
0140     if (fPups[i0] == 0)
0141       continue;
0142     // if(!fCharged[iAlgo] && fAdjust[iAlgo] && fPups[i0] > lMed) continue;
0143     if (fAdjust[iAlgo] && fPups[i0] > lMed)
0144       continue;
0145     lNRMS++;
0146     fRMS[iAlgo] += (fPups[i0] - lMed) * (fPups[i0] - lMed);
0147   }
0148   fMean[iAlgo] /= fNCount[iAlgo];
0149   if (lNRMS > 0)
0150     fRMS[iAlgo] /= lNRMS;
0151   if (fRMS[iAlgo] == 0)
0152     fRMS[iAlgo] = 1e-5;
0153   // here is the raw RMS
0154   fRMS[iAlgo] = sqrt(fRMS[iAlgo]);
0155 
0156   // some ways to do corrections to fRMS and fMedian
0157   fRMS[iAlgo] *= fRMSScaleFactor[iAlgo];
0158 
0159   if (fAdjust[iAlgo]) {
0160     //Adjust the p-value to correspond to the median
0161     int lNPV = 0;
0162     for (unsigned int i0 = 0; i0 < fPupsPV.size(); i0++)
0163       if (fPupsPV[i0] <= lMed)
0164         lNPV++;
0165     double lAdjust = double(lNPV) / double(lNPV + 0.5 * fNCount[iAlgo]);
0166     if (lAdjust > 0) {
0167       fMedian[iAlgo] -= sqrt(ROOT::Math::chisquared_quantile(lAdjust, 1.) * fRMS[iAlgo]);
0168       fRMS[iAlgo] -= sqrt(ROOT::Math::chisquared_quantile(lAdjust, 1.) * fRMS[iAlgo]);
0169     }
0170   }
0171 
0172   // fRMS_perEta[iAlgo]    *= cur_RMSEtaSF;
0173   // fMedian_perEta[iAlgo] *= cur_MedEtaSF;
0174 
0175   for (unsigned int j0 = 0; j0 < fEtaMin.size(); j0++) {
0176     fRMS_perEta[iAlgo][j0] = fRMS[iAlgo] * fRMSEtaSF[j0];
0177     fMedian_perEta[iAlgo][j0] = fMedian[iAlgo] * fMedEtaSF[j0];
0178   }
0179 }
0180 ////////////////////////////////////////////////////////////////////////////////
0181 
0182 //This code is probably a bit confusing
0183 double PuppiAlgo::compute(std::vector<double> const &iVals, double iChi2) const {
0184   if (fAlgoId[0] == -1)
0185     return 1;
0186   double lVal = 0.;
0187   double lPVal = 1.;
0188   int lNDOF = 0;
0189   for (unsigned int i0 = 0; i0 < fNAlgos; i0++) {
0190     if (fNCount[i0] == 0)
0191       return 1.;                       //in the NoPU case return 1.
0192     if (fCombId[i0] == 1 && i0 > 0) {  //Compute the previous p-value so that p-values can be multiplieed
0193       double pPVal = ROOT::Math::chisquared_cdf(lVal, lNDOF);
0194       lPVal *= pPVal;
0195       lNDOF = 0;
0196       lVal = 0;
0197     }
0198     double pVal = iVals[i0];
0199     //Special Check for any algo with log(0)
0200     if (fAlgoId[i0] == 0 && iVals[i0] == 0)
0201       pVal = cur_Med;
0202     if (fAlgoId[i0] == 3 && iVals[i0] == 0)
0203       pVal = cur_Med;
0204     if (fAlgoId[i0] == 5 && iVals[i0] == 0)
0205       pVal = cur_Med;
0206     lVal += (pVal - cur_Med) * (fabs(pVal - cur_Med)) / cur_RMS / cur_RMS;
0207     lNDOF++;
0208     if (i0 == 0 && iChi2 != 0)
0209       lNDOF++;  //Add external Chi2 to first element
0210     if (i0 == 0 && iChi2 != 0)
0211       lVal += iChi2;  //Add external Chi2 to first element
0212   }
0213   //Top it off with the last calc
0214   lPVal *= ROOT::Math::chisquared_cdf(lVal, lNDOF);
0215   return lPVal;
0216 }
0217 // ------------------------------------------------------------------------------------------
0218 void PuppiAlgo::fillDescriptionsPuppiAlgo(edm::ParameterSetDescription &desc) {
0219   edm::ParameterSetDescription puppialgos;
0220   puppialgos.add<int>("algoId", 5);
0221   puppialgos.add<bool>("useCharged", false);
0222   puppialgos.add<bool>("applyLowPUCorr", false);
0223   puppialgos.add<int>("combOpt", 5);
0224   puppialgos.add<double>("cone", .4);
0225   puppialgos.add<double>("rmsPtMin", .1);
0226   puppialgos.add<double>("rmsScaleFactor", 1.0);
0227   std::vector<edm::ParameterSet> VPSetPuppiAlgos;
0228   edm::ParameterSet puppiset;
0229   puppiset.addParameter<int>("algoId", 5);
0230   puppiset.addParameter<bool>("useCharged", false);
0231   puppiset.addParameter<bool>("applyLowPUCorr", false);
0232   puppiset.addParameter<int>("combOpt", 5);
0233   puppiset.addParameter<double>("cone", .4);
0234   puppiset.addParameter<double>("rmsPtMin", .1);
0235   puppiset.addParameter<double>("rmsScaleFactor", 1.0);
0236   VPSetPuppiAlgos.push_back(puppiset);
0237 
0238   edm::ParameterSetDescription algos;
0239   algos.addVPSet("puppiAlgos", puppialgos, VPSetPuppiAlgos);
0240   std::vector<edm::ParameterSet> VPSetAlgos;
0241   edm::ParameterSet algosset;
0242   algos.add<std::vector<double>>("etaMin", {0.});
0243   algos.add<std::vector<double>>("etaMax", {2.5});
0244   algos.add<std::vector<double>>("ptMin", {0.});
0245   algos.add<std::vector<double>>("MinNeutralPt", {0.2});
0246   algos.add<std::vector<double>>("MinNeutralPtSlope", {0.015});
0247   algos.add<std::vector<double>>("RMSEtaSF", {1.0});
0248   algos.add<std::vector<double>>("MedEtaSF", {1.0});
0249   algos.add<double>("EtaMaxExtrap", 2.0);
0250   algosset.addParameter<std::vector<double>>("etaMin", {0.});
0251   algosset.addParameter<std::vector<double>>("etaMax", {2.5});
0252   algosset.addParameter<std::vector<double>>("ptMin", {0.});
0253   algosset.addParameter<std::vector<double>>("MinNeutralPt", {0.2});
0254   algosset.addParameter<std::vector<double>>("MinNeutralPtSlope", {0.015});
0255   algosset.addParameter<std::vector<double>>("RMSEtaSF", {1.0});
0256   algosset.addParameter<std::vector<double>>("MedEtaSF", {1.0});
0257   algosset.addParameter<double>("EtaMaxExtrap", 2.0);
0258   algosset.addParameter<std::vector<edm::ParameterSet>>("puppiAlgos", VPSetPuppiAlgos);
0259   VPSetAlgos.push_back(algosset);
0260   desc.addVPSet("algos", algos, VPSetAlgos);
0261 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!