Fork me on GitHub

source: git/external/PUPPI/puppiCleanContainer.cc@ 8707eeb

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 8707eeb was fa33983, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

add PUPPI and RunPUPPI module
  • Property mode set to 100644
File size: 20.9 KB
RevLine 
[fa33983]1#include "puppiCleanContainer.hh"
2#include "Math/SpecFuncMathCore.h"
3#include "Math/ProbFunc.h"
4#include "TH2F.h"
5#include "fastjet/Selector.hh"
6
7#include <algorithm>
8
9#include "TMath.h"
10#include "Math/QuantFuncMathCore.h"
11#include "Math/SpecFuncMathCore.h"
12#include "Math/ProbFunc.h"
13
14using namespace std;
15
16// ------------- Constructor
17puppiCleanContainer::puppiCleanContainer(std::vector<RecoObj> inParticles,
18 std::vector<puppiAlgoBin> puppiAlgo,
19 float minPuppiWeight,
20 bool fUseExp){
21
22 // take the input particles
23 fRecoParticles_.clear();
24 fRecoParticles_ = inParticles;
25
26 // puppi algo
27 puppiAlgo_.clear();
28 puppiAlgo_ = puppiAlgo;
29
30 // min puppi weight
31 fMinPuppiWeight_ = minPuppiWeight;
32
33 //Clear everything
34 fPFParticles_.clear();
35 fPFchsParticles_.clear();
36 fChargedPV_.clear();
37 fChargedNoPV_.clear();
38 fPuppiWeights_.clear();
39
40 fNPV_ = 1 ;
41 fPVFrac_ = 0.;
42 fUseExp_ = fUseExp;
43
44 //Link to the RecoObjects --> loop on the input particles
45 for (unsigned int i = 0; i < fRecoParticles_.size(); i++){
46 fastjet::PseudoJet curPseudoJet;
47 curPseudoJet.reset_PtYPhiM (fRecoParticles_[i].pt,fRecoParticles_[i].eta,fRecoParticles_[i].phi,fRecoParticles_[i].m);
48 curPseudoJet.set_user_index(fRecoParticles_[i].id);
49 // fill vector of pseudojets for internal references
50 fPFParticles_.push_back(curPseudoJet);
51 if(fRecoParticles_[i].id <= 1) fPFchsParticles_.push_back(curPseudoJet); //Remove Charged particles associated to other vertex
52 if(fRecoParticles_[i].id == 1) fChargedPV_.push_back(curPseudoJet); //Take Charged particles associated to PV
53 if(fRecoParticles_[i].id == 2) fChargedNoPV_.push_back(curPseudoJet);
54 if(fRecoParticles_[i].id >= 0) fPVFrac_++ ;
55 if(fNPV_ < fRecoParticles_[i].vtxId) fNPV_ = fRecoParticles_[i].vtxId;
56
57 }
58
59 fPVFrac_ = double(fChargedPV_.size())/fPVFrac_;
60}
61
62// ------------- De-Constructor
63puppiCleanContainer::~puppiCleanContainer(){}
64
65// main function to compute puppi Event
66std::vector<fastjet::PseudoJet> puppiCleanContainer::puppiEvent(){
67
68 // output particles
69 std::vector<fastjet::PseudoJet> particles;
70 particles.clear();
71
72 std::vector<int> pPupId ;
73 std::vector<puppiParticle> partTmp ; // temp puppi particle vector; make a clone of the same particle for all the algo in which it is contained
74
75 // calculate puppi metric, RMS and mean value for all the algorithms
76 for(size_t iPuppiAlgo = 0; iPuppiAlgo < puppiAlgo_.size(); iPuppiAlgo++){
77 getRMSAvg(iPuppiAlgo,fPFParticles_,fChargedPV_); // give all the particles in the event and the charged one
78 }
79
80 int npart = 0;
81
82 // Loop on all the incoming particles
83 for(size_t iPart = 0; iPart < fPFParticles_.size(); iPart++) {
84
85 float pWeight = 1; // default weight
86 pPupId.clear();
87 pPupId = getPuppiId(fPFParticles_[iPart].pt(),fPFParticles_[iPart].eta(),puppiAlgo_); // take into account only algo eta
88
89 //////////////////////////////////////////
90 // acceptance check of the puppi algorithm
91 //////////////////////////////////////////
92
93 if(pPupId.empty()) { // out acceptance... no algorithm found
94 fPuppiWeights_.push_back(pWeight); // take the particle as it is
95 fastjet::PseudoJet curjet(pWeight*fPFParticles_[iPart].px(),pWeight*fPFParticles_[iPart].py(),pWeight*fPFParticles_[iPart].pz(),pWeight*fPFParticles_[iPart].e());
96 curjet.set_user_index(fPFParticles_[iPart].user_index());
97 particles.push_back(curjet); // fill the output collection
98 continue; //go to the next particle
99 }
100
101 ///////////////
102 // PT check //
103 ///////////////
104 for(size_t iAlgo = 0; iAlgo < pPupId.size(); iAlgo++){ // loop on all the available algo for that region
105 if(fPFParticles_.at(iPart).pt() < puppiAlgo_.at(pPupId.at(iAlgo)).fPtMin_){ // low momentum particles should be cut by puppi method
106 pWeight = 0; // if this particle is under the pT threshold of one the algorithm, put the weight as zero
107 break;
108 }
109 }
110
111 if(pWeight == 0){
112 fPuppiWeights_.push_back(0); // puppi weight is zero
113 continue;
114 }
115
116 /////////////////////////////////////
117 // fill the p-values for Z-vertex //
118 /////////////////////////////////////
119
120 double pChi2 = 0;
121 if(fUseExp_){ // use vertex-z resolution
122 //Compute an Experimental Puppi Weight with delta Z info (very simple example)
123 if(iPart <= fRecoParticles_.size() and fRecoParticles_[iPart].id == fPFParticles_.at(iPart).user_index()){
124 pChi2 = getChi2FromdZ(fRecoParticles_[iPart].dZ); // get the probability fiven the dZ of the particle wrt the leading vertex
125 if(fRecoParticles_[iPart].pfType > 3) pChi2 = 0; // not use this info for neutrals
126 }
127 }
128
129
130 /////////////////////////////////////
131 // found the particle in all the algorithm
132 /////////////////////////////////////
133 partTmp.clear();
134 for(size_t iAlgo = 0; iAlgo < pPupId.size(); iAlgo++){ // loop on all the algo found
135
136 int found = 0; // found index
137 if(fabs(fPFParticles_[iPart].user_index()) <= 1 and puppiAlgo_.at(pPupId.at(iAlgo)).fUseCharged_){ // charged or neutral from PV
138 for(size_t puppiIt = 0 ; puppiIt < puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPV_.size(); puppiIt++){ // Loop on PV particles
139 if(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPV_.at(puppiIt).fPosition_ == int(iPart)){
140 partTmp.push_back(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPV_.at(puppiIt)); // take the puppi particle
141 found = 1 ;
142 break;
143 }
144 }
145 }
146 else if ((fabs(fPFParticles_[iPart].user_index()) <= 1 and !puppiAlgo_.at(pPupId.at(iAlgo)).fUseCharged_) or fabs(fPFParticles_[iPart].user_index()) >= 2){
147 for(size_t puppiIt = 0; puppiIt < puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPU_.size(); puppiIt++){
148 if(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPU_.at(puppiIt).fPosition_ == int(iPart)){
149 partTmp.push_back(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesPU_.at(puppiIt));
150 found = 1;
151 break;
152 }
153 }
154 }
155
156 /////////////////////////////////////
157 // means that is inside the NULL vector for some reasons
158 /////////////////////////////////////
159
160 if(found == 0){
161 for(size_t puppiIt = 0; puppiIt < puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesNULL_.size(); puppiIt++){
162 if(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesNULL_.at(puppiIt).fPosition_ == int(iPart)){
163 partTmp.push_back(puppiAlgo_.at(pPupId.at(iAlgo)).fPuppiParticlesNULL_.at(puppiIt));
164 found = 1 ;
165 break;
166 }
167 }
168 }
169 }
170
171 if(partTmp.size() != pPupId.size()){ // not found the particle in one of the algorithms
172 pWeight = 1 ;
173 fPuppiWeights_.push_back(pWeight);
174 fastjet::PseudoJet curjet( pWeight*fPFParticles_[iPart].px(), pWeight*fPFParticles_[iPart].py(), pWeight*fPFParticles_[iPart].pz(), pWeight*fPFParticles_[iPart].e());
175 curjet.set_user_index(fPFParticles_[iPart].user_index());
176 particles.push_back(curjet); // by default is one, so 4V is not chaged
177 continue;
178 }
179
180 /////////////////////////////////////
181 //Check the Pval
182 /////////////////////////////////////
183
184 bool badPVal = false ;
185 for(size_t iPuppi = 0; iPuppi < partTmp.size() ; iPuppi++){
186 if(partTmp.at(iPuppi).fPval_ == -999){ // if the default is found as PVal, leave the particle as it is in the output
187 pWeight = 1 ;
188 badPVal = true ;
189 }
190 }
191
192 if(badPVal){
193 fPuppiWeights_.push_back(pWeight);
194 fastjet::PseudoJet curjet( pWeight*fPFParticles_[iPart].px(), pWeight*fPFParticles_[iPart].py(), pWeight*fPFParticles_[iPart].pz(), pWeight*fPFParticles_[iPart].e());
195 curjet.set_user_index(fPFParticles_[iPart].user_index());
196 particles.push_back(curjet);
197 continue;
198
199 }
200
201
202 // compute combining the weight for all the algorithm
203 pWeight = compute(pChi2,partTmp,puppiAlgo_,pPupId);
204
205 //Basic Weight Checks
206 if( std::isinf(pWeight) || std::isnan(pWeight)){
207 std::cerr << "====> Weight is nan : pt " << fPFParticles_[iPart].pt() << " -- eta : " << fPFParticles_[iPart].eta() << " -- id : " << fPFParticles_[iPart].user_index() << std::endl;
208 pWeight = 1; // set the default to avoid problems
209 }
210
211 //Basic Cuts
212 if(pWeight < fMinPuppiWeight_) pWeight = 0; //==> Elminate the low Weight stuff
213
214 //threshold cut on the neutral Pt
215 for(size_t iPuppi = 0; iPuppi < pPupId.size(); iPuppi++){
216 if(fPFParticles_[iPart].user_index() == 1 && puppiAlgo_.at(pPupId.at(iPuppi)).fApplyCHS_ ) pWeight = 1; // charged from LV
217 if(fPFParticles_[iPart].user_index() == 2 && puppiAlgo_.at(pPupId.at(iPuppi)).fApplyCHS_ ) pWeight = 0; // charged from PU
218 if(pWeight*fPFParticles_[iPart].pt() < getNeutralPtCut(puppiAlgo_.at(pPupId.at(iPuppi)).fNeutralMinE_,puppiAlgo_.at(pPupId.at(iPuppi)).fNeutralPtSlope_,fNPV_) && fPFParticles_[iPart].user_index() == 0 ) // if don't pass one of the algo neutral pt condition the particle is cut
219 pWeight = 0;
220 }
221
222 fPuppiWeights_.push_back(pWeight); // push back the weight
223
224 //Now get rid of the thrown out weights for the particle collection
225 if(pWeight == 0) continue; // if zero don't fill the particle in the output
226 npart++;
227
228 //Produce
229 fastjet::PseudoJet curjet( pWeight*fPFParticles_[iPart].px(), pWeight*fPFParticles_[iPart].py(), pWeight*fPFParticles_[iPart].pz(), pWeight*fPFParticles_[iPart].e());
230 curjet.set_user_index(iPart);
231 particles.push_back(curjet);
232
233 }
234
235 return particles;
236
237}
238
239// compute puppi metric, RMS and median for PU particle for each algo
240void puppiCleanContainer::getRMSAvg(const int & iPuppiAlgo, std::vector<fastjet::PseudoJet> & particlesAll, std::vector<fastjet::PseudoJet> &chargedPV) {
241
242 std::vector<puppiParticle> puppiParticles; // puppi particles to be set for a specific algo
243 puppiParticles.clear();
244
245 // Loop on all the particles of the event
246
247 for(size_t iPart = 0; iPart < particlesAll.size(); iPart++ ) {
248
249 float pVal = -999;
250 bool pPupId = isGoodPuppiId(particlesAll[iPart].pt(),particlesAll[iPart].eta(),puppiAlgo_.at(iPuppiAlgo)); // get the puppi id algo asaf of eta and phi of the particle
251 // does not exsist and algorithm for this particle, store -999 as pVal
252 if(pPupId == false) continue;
253 // apply CHS in puppi metric computation -> use only LV hadrons to compute the metric for each particle
254 if(puppiAlgo_.at(iPuppiAlgo).fUseCharged_)
255 pVal = goodVar(particlesAll[iPart], chargedPV, puppiAlgo_.at(iPuppiAlgo).fMetricId_,puppiAlgo_.at(iPuppiAlgo).fConeSize_);
256 else if(!puppiAlgo_.at(iPuppiAlgo).fUseCharged_)
257 pVal = goodVar(particlesAll[iPart], particlesAll, puppiAlgo_.at(iPuppiAlgo).fMetricId_,puppiAlgo_.at(iPuppiAlgo).fConeSize_);
258
259 // fill the value
260 if(std::isnan(pVal) || std::isinf(pVal)) std::cout << "====> Value is Nan " << pVal << " == " << particlesAll[iPart].pt() << " -- " << particlesAll[iPart].eta() << std::endl;
261 if(std::isnan(pVal) || std::isinf(pVal)) continue;
262
263 puppiParticles.push_back(puppiParticle(particlesAll.at(iPart).pt(),particlesAll.at(iPart).eta(),pVal,particlesAll.at(iPart).user_index(),iPart));
264 }
265
266 // set the puppi particles for the algorithm
267 puppiAlgo_.at(iPuppiAlgo).setPuppiParticles(puppiParticles);
268 // compute RMS, median and mean value
269 computeMedRMS(iPuppiAlgo);
270
271}
272
273
274float puppiCleanContainer::goodVar(const fastjet::PseudoJet & particle, const std::vector<fastjet::PseudoJet> & particleAll, const int & pPupId, const float & coneSize) {
275 float lPup = 0;
276 lPup = var_within_R(pPupId,particleAll,particle,coneSize);
277 return lPup;
278}
279
280float puppiCleanContainer::var_within_R(const int & pPupId, const vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet& centre, const float & R){
281
282 if(pPupId == -1) return 1;
283 fastjet::Selector sel = fastjet::SelectorCircle(R);
284 sel.set_reference(centre);
285 std::vector<fastjet::PseudoJet> near_particles = sel(particles);
286 float var = 0;
287
288 for(size_t iPart = 0; iPart < near_particles.size(); iPart++){
289
290 double pDEta = near_particles[iPart].eta()-centre.eta();
291 double pDPhi = fabs(near_particles[iPart].phi()-centre.phi());
292 if(pDPhi > 2.*3.14159265-pDPhi) pDPhi = 2.*3.14159265-pDPhi;
293 double pDR = sqrt(pDEta*pDEta+pDPhi*pDPhi);
294
295 if(pDR < 0.0001) continue;
296 if(pDR == 0) continue;
297
298 if(pPupId == 0) var += (near_particles[iPart].pt()/(pDR*pDR));
299 if(pPupId == 1) var += near_particles[iPart].pt();
300 if(pPupId == 2) var += (1./pDR)*(1./pDR);
301 if(pPupId == 3) var += (1./pDR)*(1./pDR);
302 if(pPupId == 4) var += near_particles[iPart].pt();
303 if(pPupId == 5) var += (near_particles[iPart].pt()/pDR)*(near_particles[iPart].pt()/pDR);
304 }
305
306 if(pPupId == 0 && var != 0) var = log(var);
307 if(pPupId == 3 && var != 0) var = log(var);
308 if(pPupId == 5 && var != 0) var = log(var);
309 return var;
310
311}
312
313
314float puppiCleanContainer::pt_within_R(const std::vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet & centre, const float & R){
315
316 fastjet::Selector sel = fastjet::SelectorCircle(R);
317 sel.set_reference(centre);
318 std::vector<fastjet::PseudoJet> near_particles = sel(particles);
319 double answer = 0.0;
320 for(size_t iPart = 0; iPart<near_particles.size(); iPart++){
321 answer += near_particles[iPart].pt();
322 }
323
324 return answer;
325}
326
327fastjet::PseudoJet puppiCleanContainer::flow_within_R(const vector<fastjet::PseudoJet> & particles, const fastjet::PseudoJet& centre, const float & R){
328
329 fastjet::Selector sel = fastjet::SelectorCircle(R);
330 sel.set_reference(centre);
331 std::vector<fastjet::PseudoJet> near_particles = sel(particles);
332 fastjet::PseudoJet flow;
333 for(unsigned int i=0; i<near_particles.size(); i++){
334 flow += near_particles[i];
335 }
336 return flow;
337
338}
339
340// compute median, mean value and RMS for puppi
341void puppiCleanContainer::computeMedRMS(const int & puppiAlgo) {
342
343 if(puppiAlgo > int(puppiAlgo_.size()) ) return;
344 if(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size() == 0) return;
345
346 // sort in pVal increasing order
347 std::sort(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.begin(),puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.end(),puppiValSort());
348
349 // if apply correction
350 float lCorr = 1.;
351 if(puppiAlgo_.at(puppiAlgo).fApplyLowPUCorr_) lCorr *= 1.-fPVFrac_;
352
353 // count the position of the last particle with pval zero coming from PU
354 int lNum0 = 0;
355 for(size_t i0 = 0; i0 < puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size(); i0++) {
356 if(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_[i0].fPval_ == 0) lNum0 = i0;
357 }
358
359 // take the median value on PU particles
360 int lNHalfway = lNum0 + int(float(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size()-lNum0)*0.50*lCorr);
361 puppiAlgo_.at(puppiAlgo).fMedian_ = puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(lNHalfway).fPval_;
362 float lMed = puppiAlgo_.at(puppiAlgo).fMedian_; //Just to make the readability easier
363
364 // take the RMS
365 int lNRMS = 0;
366 for(size_t i0 = 0; i0 < puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size(); i0++) {
367 puppiAlgo_.at(puppiAlgo).fMean_ += puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_;
368 if(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_ == 0) continue;
369 if(!puppiAlgo_.at(puppiAlgo).fUseCharged_ && puppiAlgo_.at(puppiAlgo).fApplyLowPUCorr_ && puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_ > lMed) continue;
370 lNRMS++;
371 puppiAlgo_.at(puppiAlgo).fRMS_ += (puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_-lMed)*( puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.at(i0).fPval_-lMed);
372 }
373
374 puppiAlgo_.at(puppiAlgo).fMean_ /= puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size();
375 if(lNRMS > 0) puppiAlgo_.at(puppiAlgo).fRMS_/=lNRMS;
376 if(puppiAlgo_.at(puppiAlgo).fRMS_ == 0) puppiAlgo_.at(puppiAlgo).fRMS_ = 1e-5;
377 puppiAlgo_.at(puppiAlgo).fRMS_ = sqrt(puppiAlgo_.at(puppiAlgo).fRMS_);
378 puppiAlgo_.at(puppiAlgo).fRMS_ *= puppiAlgo_.at(puppiAlgo).fRMSScaleFactor_;
379
380 if(!puppiAlgo_.at(puppiAlgo).fApplyLowPUCorr_) return;
381
382 //Adjust the p-value to correspond to the median
383 std::sort(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_.begin(),puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_.end(),puppiValSort());
384 int lNPV = 0;
385 for(size_t i0 = 0; i0 < puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_.size(); i0++){
386 if(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_[i0].fPval_ <= lMed ) lNPV++;
387 }
388
389 // it helps in puppi the median value close to the mean one when a lot of pval 0 are present
390 float lAdjust = 1.5*float(lNPV)/float(puppiAlgo_.at(puppiAlgo).fPuppiParticlesPV_.size()+puppiAlgo_.at(puppiAlgo).fPuppiParticlesPU_.size());
391 if(lAdjust > 0) puppiAlgo_.at(puppiAlgo).fMedian_ -= sqrt(ROOT::Math::chisquared_quantile(lAdjust,1.)*puppiAlgo_.at(puppiAlgo).fRMS_);
392
393}
394
395float puppiCleanContainer::getNeutralPtCut(const float & fNeutralMinE, const float & fNeutralPtSlope, const int & fNPV) {
396 return fNeutralMinE + fNPV * fNeutralPtSlope;
397}
398
399// take the type of algorithm : return a vector since more than one algo can be defined for the same eta region
400std::vector<int> puppiCleanContainer::getPuppiId(const float & pt, const float & eta, const std::vector<puppiAlgoBin> & puppiAlgos){
401 std::vector<int> PuppiId ;
402 for(size_t iPuppiAlgo = 0; iPuppiAlgo < puppiAlgos.size() ; iPuppiAlgo++){
403 if(fabs(eta) <= puppiAlgos[iPuppiAlgo].fEtaMin_) continue;
404 if(fabs(eta) > puppiAlgos[iPuppiAlgo].fEtaMax_) continue;
405 PuppiId.push_back(int(iPuppiAlgo));
406 }
407 return PuppiId;
408}
409
410//check if a particle is good for an Algo definition
411bool puppiCleanContainer::isGoodPuppiId(const float & pt, const float & eta, const puppiAlgoBin & puppiAlgo){
412 if(fabs(eta) <= puppiAlgo.fEtaMin_) return false;
413 if(fabs(eta) > puppiAlgo.fEtaMax_) return false;
414 return true;
415
416}
417
418
419
420// ----------------------
421float puppiCleanContainer::compute(const float & chi2, const std::vector<puppiParticle> & particles, const std::vector<puppiAlgoBin> & puppiAlgos, const std::vector<int> & pPupId) {
422
423 if(particles.size() != pPupId.size() ) return 0; // default check
424
425 float lVal = 0.;
426 float lPVal = 1.;
427 int lNDOF = 0;
428
429 for( size_t iAlgo = 0; iAlgo < pPupId.size(); iAlgo++){
430
431 if(puppiAlgos.at(pPupId.at(iAlgo)).fPuppiParticlesPU_.size() + puppiAlgos.at(pPupId.at(iAlgo)).fPuppiParticlesPV_.size() == 0) return 1;
432
433 if(iAlgo > 0 ){
434 float pPVal = ROOT::Math::chisquared_cdf(lVal,lNDOF); // take a chi2 value since the should be multiplied (multiply weight and not summing chi2)
435 lPVal *= pPVal;
436 lNDOF = 0;
437 lVal = 0;
438 }
439
440 if(puppiAlgos.at(pPupId.at(iAlgo)).fMetricId_ == -1) continue;
441
442 float pVal = particles.at(iAlgo).fPval_ ;
443
444 if(puppiAlgos.at(pPupId.at(iAlgo)).fMetricId_ == 0 && pVal == 0) pVal = puppiAlgos.at(pPupId.at(iAlgo)).fMedian_;
445 if(puppiAlgos.at(pPupId.at(iAlgo)).fMetricId_ == 3 && pVal == 0) pVal = puppiAlgos.at(pPupId.at(iAlgo)).fMedian_;
446 if(puppiAlgos.at(pPupId.at(iAlgo)).fMetricId_ == 5 && pVal == 0) pVal = puppiAlgos.at(pPupId.at(iAlgo)).fMedian_;
447
448 lVal += (pVal-puppiAlgos.at(pPupId.at(iAlgo)).fMedian_)*(fabs(pVal-puppiAlgos.at(pPupId.at(iAlgo)).fMedian_))/puppiAlgos.at(pPupId.at(iAlgo)).fRMS_/puppiAlgos.at(pPupId.at(iAlgo)).fRMS_;
449 lNDOF++;
450 if(chi2 != 0) lNDOF++;
451 if(chi2 != 0) lVal+=chi2; //Add external Chi2 to first element
452 }
453
454 lPVal *= ROOT::Math::chisquared_cdf(lVal,lNDOF);
455 return lPVal;
456
457}
458
459float puppiCleanContainer::getChi2FromdZ(float iDZ) {
460 //We need to obtain prob of PU + (1-Prob of LV)
461 // Prob(LV) = Gaus(dZ,sigma) where sigma = 1.5mm (its really more like 1mm)
462 //double lProbLV = ROOT::Math::normal_cdf_c(fabs(iDZ),0.2)*2.; //*2 is to do it double sided
463 //Take iDZ to be corrected by sigma already
464 double lProbLV = ROOT::Math::normal_cdf_c(fabs(iDZ),1.)*2.; //*2 is to do it double sided
465 double lProbPU = 1-lProbLV;
466 if(lProbPU <= 0) lProbPU = 1e-16; //Quick Trick to through out infs
467 if(lProbPU >= 0) lProbPU = 1-1e-16; //Ditto
468 double lChi2PU = TMath::ChisquareQuantile(lProbPU,1);
469 lChi2PU*=lChi2PU;
470 return lChi2PU;
471}
Note: See TracBrowser for help on using the repository browser.