| 1 | #include "PuppiContainer.hh"
|
|---|
| 2 | #include "external/fastjet/internal/base.hh"
|
|---|
| 3 | #include "external/fastjet/Selector.hh"
|
|---|
| 4 | //#include "external/fastjet/FunctionOfPseudoJet.hh"
|
|---|
| 5 | #include "Math/ProbFunc.h"
|
|---|
| 6 | #include "TMath.h"
|
|---|
| 7 | #include <iostream>
|
|---|
| 8 | #include <math.h>
|
|---|
| 9 |
|
|---|
| 10 | PuppiContainer::PuppiContainer(bool iApplyCHS, bool iUseExp,double iPuppiWeightCut,std::vector<AlgoObj> &iAlgos) {
|
|---|
| 11 | fApplyCHS = iApplyCHS;
|
|---|
| 12 | fUseExp = iUseExp;
|
|---|
| 13 | fPuppiWeightCut = iPuppiWeightCut;
|
|---|
| 14 | fNAlgos = iAlgos.size();
|
|---|
| 15 | for(unsigned int i0 = 0; i0 < iAlgos.size(); i0++) {
|
|---|
| 16 | PuppiAlgo pPuppiConfig(iAlgos[i0]);
|
|---|
| 17 | fPuppiAlgo.push_back(pPuppiConfig);
|
|---|
| 18 | }
|
|---|
| 19 | }
|
|---|
| 20 |
|
|---|
| 21 | void PuppiContainer::initialize(const std::vector<RecoObj> &iRecoObjects) {
|
|---|
| 22 | //Clear everything
|
|---|
| 23 | fRecoParticles.resize(0);
|
|---|
| 24 | fPFParticles .resize(0);
|
|---|
| 25 | fChargedPV .resize(0);
|
|---|
| 26 | fPupParticles .resize(0);
|
|---|
| 27 | fWeights .resize(0);
|
|---|
| 28 | fVals.resize(0);
|
|---|
| 29 | //fChargedNoPV.resize(0);
|
|---|
| 30 | //Link to the RecoObjects
|
|---|
| 31 | fPVFrac = 0.;
|
|---|
| 32 | fNPV = 1.;
|
|---|
| 33 | fRecoParticles = iRecoObjects;
|
|---|
| 34 | for (unsigned int i = 0; i < fRecoParticles.size(); i++){
|
|---|
| 35 | fastjet::PseudoJet curPseudoJet;
|
|---|
| 36 | curPseudoJet.reset_PtYPhiM(fRecoParticles[i].pt,fRecoParticles[i].eta,fRecoParticles[i].phi,fRecoParticles[i].m);
|
|---|
| 37 | if(fRecoParticles[i].id == 0 or fRecoParticles[i].charge == 0) curPseudoJet.set_user_index(0); // zero is neutral hadron
|
|---|
| 38 | if(fRecoParticles[i].id == 1 and fRecoParticles[i].charge != 0) curPseudoJet.set_user_index(fRecoParticles[i].charge); // from PV use the
|
|---|
| 39 | if(fRecoParticles[i].id == 2 and fRecoParticles[i].charge != 0) curPseudoJet.set_user_index(fRecoParticles[i].charge+5); // from NPV use the charge as key +5 as key // fill vector of pseudojets for internal references
|
|---|
| 40 | fPFParticles.push_back(curPseudoJet);
|
|---|
| 41 | //Take Charged particles associated to PV
|
|---|
| 42 | if(fabs(fRecoParticles[i].id) == 1) fChargedPV.push_back(curPseudoJet);
|
|---|
| 43 | if(fabs(fRecoParticles[i].id) >= 1 ) fPVFrac+=1.;
|
|---|
| 44 | if(fNPV < fRecoParticles[i].vtxId) fNPV = fRecoParticles[i].vtxId;
|
|---|
| 45 | }
|
|---|
| 46 | //if(fNPV < 10) fNPV = 80.;
|
|---|
| 47 | if(fPVFrac != 0) { fPVFrac = double(fChargedPV.size())/fPVFrac;}
|
|---|
| 48 | else { fPVFrac = 0;}
|
|---|
| 49 | }
|
|---|
| 50 | PuppiContainer::~PuppiContainer(){}
|
|---|
| 51 |
|
|---|
| 52 | double PuppiContainer::goodVar(fastjet::PseudoJet &iPart,std::vector<fastjet::PseudoJet> &iParts, int iOpt,double iRCone) {
|
|---|
| 53 | double lPup = 0;
|
|---|
| 54 | lPup = var_within_R(iOpt,iParts,iPart,iRCone);
|
|---|
| 55 | return lPup;
|
|---|
| 56 | }
|
|---|
| 57 | double PuppiContainer::var_within_R(int iId, const vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet& centre, double R){
|
|---|
| 58 | if(iId == -1) return 1;
|
|---|
| 59 | fastjet::Selector sel = fastjet::SelectorCircle(R);
|
|---|
| 60 | sel.set_reference(centre);
|
|---|
| 61 | vector<fastjet::PseudoJet> near_particles = sel(particles);
|
|---|
| 62 | double var = 0;
|
|---|
| 63 | //double lSumPt = 0;
|
|---|
| 64 | //if(iId == 1) for(unsigned int i=0; i<near_particles.size(); i++) lSumPt += near_particles[i].pt();
|
|---|
| 65 | for(unsigned int i=0; i<near_particles.size(); i++){
|
|---|
| 66 | double pDEta = near_particles[i].eta()-centre.eta();
|
|---|
| 67 | double pDPhi = fabs(near_particles[i].phi()-centre.phi());
|
|---|
| 68 | if(pDPhi > 2.*3.14159265-pDPhi) pDPhi = 2.*3.14159265-pDPhi;
|
|---|
| 69 | double pDR2 = pDEta*pDEta+pDPhi*pDPhi;
|
|---|
| 70 | if(std::abs(pDR2) < 0.0001) continue;
|
|---|
| 71 | if(iId == 0) var += (near_particles[i].pt()/pDR2);
|
|---|
| 72 | if(iId == 1) var += near_particles[i].pt();
|
|---|
| 73 | if(iId == 2) var += (1./pDR2);
|
|---|
| 74 | if(iId == 3) var += (1./pDR2);
|
|---|
| 75 | if(iId == 4) var += near_particles[i].pt();
|
|---|
| 76 | if(iId == 5) var += (near_particles[i].pt()*(near_particles[i].pt()/pDR2));
|
|---|
| 77 | }
|
|---|
| 78 | if(iId == 1) var += centre.pt(); //Sum in a cone
|
|---|
| 79 | if(iId == 0 && var != 0) var = log(var);
|
|---|
| 80 | if(iId == 3 && var != 0) var = log(var);
|
|---|
| 81 | if(iId == 5 && var != 0) var = log(var);
|
|---|
| 82 | return var;
|
|---|
| 83 | }
|
|---|
| 84 | //In fact takes the median not the average
|
|---|
| 85 | void PuppiContainer::getRMSAvg(int iOpt,std::vector<fastjet::PseudoJet> &iConstits,std::vector<fastjet::PseudoJet> &iParticles,std::vector<fastjet::PseudoJet> &iChargedParticles) {
|
|---|
| 86 | for(unsigned int i0 = 0; i0 < iConstits.size(); i0++ ) {
|
|---|
| 87 | double pVal = -1;
|
|---|
| 88 | //Calculate the Puppi Algo to use
|
|---|
| 89 | int pPupId = getPuppiId(iConstits[i0].pt(),iConstits[i0].eta());
|
|---|
| 90 | if(fPuppiAlgo[pPupId].numAlgos() <= iOpt) pPupId = -1;
|
|---|
| 91 | if(pPupId == -1) {fVals.push_back(-1); continue;}
|
|---|
| 92 | //Get the Puppi Sub Algo (given iteration)
|
|---|
| 93 | int pAlgo = fPuppiAlgo[pPupId].algoId (iOpt);
|
|---|
| 94 | bool pCharged = fPuppiAlgo[pPupId].isCharged(iOpt);
|
|---|
| 95 | double pCone = fPuppiAlgo[pPupId].coneSize (iOpt);
|
|---|
| 96 | //Compute the Puppi Metric
|
|---|
| 97 | if(!pCharged) pVal = goodVar(iConstits[i0],iParticles ,pAlgo,pCone);
|
|---|
| 98 | if( pCharged) pVal = goodVar(iConstits[i0],iChargedParticles,pAlgo,pCone);
|
|---|
| 99 | fVals.push_back(pVal);
|
|---|
| 100 | if(std::isnan(pVal) || std::isinf(pVal)) cerr << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt() << " -- " << iConstits[i0].eta() << endl;
|
|---|
| 101 | if(std::isnan(pVal) || std::isinf(pVal)) continue;
|
|---|
| 102 | fPuppiAlgo[pPupId].add(iConstits[i0],pVal,iOpt);
|
|---|
| 103 | /* eta extrapolation code
|
|---|
| 104 | for(int i1 = 0; i1 < fNAlgos; i1++){
|
|---|
| 105 | pAlgo = fPuppiAlgo[i1].algoId (iOpt);
|
|---|
| 106 | pCharged = fPuppiAlgo[i1].isCharged(iOpt);
|
|---|
| 107 | pCone = fPuppiAlgo[i1].coneSize (iOpt);
|
|---|
| 108 | double curVal = -1;
|
|---|
| 109 | if(!pCharged) curVal = goodVar(iConstits[i0],iParticles ,pAlgo,pCone);
|
|---|
| 110 | if( pCharged) curVal = goodVar(iConstits[i0],iChargedParticles,pAlgo,pCone);
|
|---|
| 111 | //std::cout << "i1 = " << i1 << ", curVal = " << curVal << ", eta = " << iConstits[i0].eta() << ", pupID = " << pPupId << std::endl;
|
|---|
| 112 | fPuppiAlgo[i1].add(iConstits[i0],curVal,iOpt);
|
|---|
| 113 | }
|
|---|
| 114 | */
|
|---|
| 115 | }
|
|---|
| 116 | for(int i0 = 0; i0 < fNAlgos; i0++) fPuppiAlgo[i0].computeMedRMS(iOpt,fPVFrac);
|
|---|
| 117 | }
|
|---|
| 118 | int PuppiContainer::getPuppiId(const float &iPt,const float &iEta) {
|
|---|
| 119 | int lId = -1;
|
|---|
| 120 | for(int i0 = 0; i0 < fNAlgos; i0++) {
|
|---|
| 121 | if(fabs(iEta) < fPuppiAlgo[i0].etaMin()) continue;
|
|---|
| 122 | if(fabs(iEta) > fPuppiAlgo[i0].etaMax()) continue;
|
|---|
| 123 | if(iPt < fPuppiAlgo[i0].ptMin()) continue;
|
|---|
| 124 | lId = i0;
|
|---|
| 125 | break;
|
|---|
| 126 | }
|
|---|
| 127 | return lId;
|
|---|
| 128 | }
|
|---|
| 129 | double PuppiContainer::getChi2FromdZ(double iDZ) {
|
|---|
| 130 | //We need to obtain prob of PU + (1-Prob of LV)
|
|---|
| 131 | // Prob(LV) = Gaus(dZ,sigma) where sigma = 1.5mm (its really more like 1mm)
|
|---|
| 132 | //double lProbLV = ROOT::Math::normal_cdf_c(fabs(iDZ),0.2)*2.; //*2 is to do it double sided
|
|---|
| 133 | //Take iDZ to be corrected by sigma already
|
|---|
| 134 | double lProbLV = ROOT::Math::normal_cdf_c(fabs(iDZ),1.)*2.; //*2 is to do it double sided
|
|---|
| 135 | double lProbPU = 1-lProbLV;
|
|---|
| 136 | if(lProbPU <= 0) lProbPU = 1e-16; //Quick Trick to through out infs
|
|---|
| 137 | if(lProbPU >= 0) lProbPU = 1-1e-16; //Ditto
|
|---|
| 138 | double lChi2PU = TMath::ChisquareQuantile(lProbPU,1);
|
|---|
| 139 | lChi2PU*=lChi2PU;
|
|---|
| 140 | return lChi2PU;
|
|---|
| 141 | }
|
|---|
| 142 | const std::vector<double> PuppiContainer::puppiWeights() {
|
|---|
| 143 | fPupParticles .resize(0);
|
|---|
| 144 | fWeights .resize(0);
|
|---|
| 145 | fVals .resize(0);
|
|---|
| 146 | for(int i0 = 0; i0 < fNAlgos; i0++) fPuppiAlgo[i0].reset();
|
|---|
| 147 |
|
|---|
| 148 | int lNMaxAlgo = 1;
|
|---|
| 149 | for(int i0 = 0; i0 < fNAlgos; i0++) lNMaxAlgo = TMath::Max(fPuppiAlgo[i0].numAlgos(),lNMaxAlgo);
|
|---|
| 150 | //Run through all compute mean and RMS
|
|---|
| 151 | int lNParticles = fRecoParticles.size();
|
|---|
| 152 | for(int i0 = 0; i0 < lNMaxAlgo; i0++) {
|
|---|
| 153 | getRMSAvg(i0,fPFParticles,fPFParticles,fChargedPV);
|
|---|
| 154 | }
|
|---|
| 155 | std::vector<double> pVals;
|
|---|
| 156 | for(int i0 = 0; i0 < lNParticles; i0++) {
|
|---|
| 157 | //Refresh
|
|---|
| 158 | pVals.clear();
|
|---|
| 159 | double pWeight = 1;
|
|---|
| 160 | //Get the Puppi Id and if ill defined move on
|
|---|
| 161 | int pPupId = getPuppiId(fRecoParticles[i0].pt,fRecoParticles[i0].eta);
|
|---|
| 162 | if(pPupId == -1) {
|
|---|
| 163 | fWeights .push_back(pWeight);
|
|---|
| 164 | continue;
|
|---|
| 165 | }
|
|---|
| 166 | // fill the p-values
|
|---|
| 167 | double pChi2 = 0;
|
|---|
| 168 | if(fUseExp){
|
|---|
| 169 | //Compute an Experimental Puppi Weight with delta Z info (very simple example)
|
|---|
| 170 | pChi2 = getChi2FromdZ(fRecoParticles[i0].dZ);
|
|---|
| 171 | //Now make sure Neutrals are not set
|
|---|
| 172 | if(fRecoParticles[i0].pfType > 3) pChi2 = 0;
|
|---|
| 173 | }
|
|---|
| 174 | //Fill and compute the PuppiWeight
|
|---|
| 175 | int lNAlgos = fPuppiAlgo[pPupId].numAlgos();
|
|---|
| 176 | for(int i1 = 0; i1 < lNAlgos; i1++) pVals.push_back(fVals[lNParticles*i1+i0]);
|
|---|
| 177 | pWeight = fPuppiAlgo[pPupId].compute(pVals,pChi2);
|
|---|
| 178 | //Apply the CHS weights
|
|---|
| 179 | if(fRecoParticles[i0].id == 1 && fApplyCHS ) pWeight = 1;
|
|---|
| 180 | if(fRecoParticles[i0].id == 2 && fApplyCHS ) pWeight = 0;
|
|---|
| 181 | //Basic Weight Checks
|
|---|
| 182 | if(std::isnan(pWeight)) std::cerr << "====> Weight is nan : pt " << fRecoParticles[i0].pt << " -- eta : " << fRecoParticles[i0].eta << " -- Value" << fVals[i0] << " -- id : " << fRecoParticles[i0].id << " -- NAlgos: " << lNAlgos << std::endl;
|
|---|
| 183 | //Basic Cuts
|
|---|
| 184 | if(pWeight < fPuppiWeightCut) pWeight = 0; //==> Elminate the low Weight stuff
|
|---|
| 185 | if(pWeight*fPFParticles[i0].pt() < fPuppiAlgo[pPupId].neutralPt(fNPV) && fRecoParticles[i0].id == 0 ) pWeight = 0; //threshold cut on the neutral Pt
|
|---|
| 186 | fWeights .push_back(pWeight);
|
|---|
| 187 | //Now get rid of the thrown out weights for the particle collection
|
|---|
| 188 | if(pWeight == 0) continue;
|
|---|
| 189 | //Produce
|
|---|
| 190 | fastjet::PseudoJet curjet( pWeight*fPFParticles[i0].px(), pWeight*fPFParticles[i0].py(), pWeight*fPFParticles[i0].pz(), pWeight*fPFParticles[i0].e());
|
|---|
| 191 | curjet.set_user_index(i0);//fRecoParticles[i0].id);
|
|---|
| 192 | fPupParticles.push_back(curjet);
|
|---|
| 193 | }
|
|---|
| 194 | return fWeights;
|
|---|
| 195 | }
|
|---|
| 196 |
|
|---|
| 197 |
|
|---|
