Fork me on GitHub

source: git/external/PUPPI/PuppiContainer.cc@ 0e6b14e

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

fix PUPPI includes

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