ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitPhysics/Utils/src/JetTools.cc

Comparing UserCode/MitPhysics/Utils/src/JetTools.cc (file contents):
Revision 1.3 by ceballos, Fri Jul 30 20:23:31 2010 UTC vs.
Revision 1.33 by pharris, Thu Mar 21 17:50:13 2013 UTC

#	Line 1 | Line 1
1		#include "MitPhysics/Utils/interface/JetTools.h"
2	+	#include <algorithm>
3	+	#include <vector>
4	+	#include "TMatrixDSym.h"
5	+	#include "TMatrixDSymEigen.h"
6
7		ClassImp(mithep::JetTools)
8
9		using namespace mithep;
10	+
11
12		JetTools::JetTools()
13		{
#	Line 15 | Line 20 | JetTools::~JetTools()
20		}
21
22		//Remember to remove the signal from particles before inputting into the function
23	<	Double_t JetTools::NJettiness(const ParticleOArr *particles, const JetOArr *jets, bool UseQ, Double_t Y){
23	>	Double_t JetTools::NJettiness(const ParticleOArr *particles, const JetOArr *jets, double Q, double Y){
24		if(particles->GetEntries() <= 0) return 0.0;
25
26		Double_t fval = 0.0;
#	Line 27 | Line 32 | Double_t JetTools::NJettiness(const Part
32		for(int j=0;j<int(jets->GetEntries());j++){
33		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
34		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-particles->At(i)->Eta()))
35	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi()))));
35	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi())))));
36		}
37		fval = fval + fvalpart;
38		}
39
40	<	if(UseQ == kTRUE) fval = fval / particles->At(0)->Pt();
40	>	fval = fval / Q;
41
42		return fval;
43		}
44
45	<	Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, bool UseQ, Double_t Y){
45	>	Double_t JetTools::NJettiness(const PFCandidateOArr *pfCandidates, const JetOArr *jets, double Q, double Y){
46	>	if(pfCandidates->GetEntries() <= 0) return 0.0;
47	>
48	>	Double_t fval = 0.0;
49	>	Double_t fvalpart;
50	>
51	>	for(int i=0;i<int(pfCandidates->GetEntries());i++){
52	>	fvalpart = (pfCandidates->At(i)->Pt()) * TMath::Exp(-TMath::Abs(pfCandidates->At(i)->Eta()-Y));
53	>
54	>	for(int j=0;j<int(jets->GetEntries());j++){
55	>	fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
56	>	(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-pfCandidates->At(i)->Eta()))
57	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi())))));
58	>	}
59	>	fval = fval + fvalpart;
60	>	}
61	>
62	>	fval = fval / Q;
63	>
64	>	return fval;
65	>	}
66	>
67	>	Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, double Q, double Y){
68		if(tracks->GetEntries() <= 0) return 0.0;
69
70		Double_t fval = 0.0;
#	Line 49 | Line 76 | Double_t JetTools::NJettiness(const Trac
76		for(int j=0;j<int(jets->GetEntries());j++){
77		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
78		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-tracks->At(i)->Eta()))
79	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi()))));
79	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi())))));
80		}
81		fval = fval + fvalpart;
82		}
83
84	<	if(UseQ == kTRUE) fval = fval / tracks->At(0)->Pt();
84	>	fval = fval / Q;
85
86		return fval;
87		}
88
89	<	Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, bool UseQ, Double_t Y){
89	>	Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, double Q, double Y){
90		if(jetsS->GetEntries() <= 0) return 0.0;
91
92		Double_t fval = 0.0;
#	Line 71 | Line 98 | Double_t JetTools::NJettiness(const JetO
98		for(int j=0;j<int(jets->GetEntries());j++){
99		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
100		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-jetsS->At(i)->Eta()))
101	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi()))));
101	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi())))));
102		}
103		fval = fval + fvalpart;
104		}
105
106	<	if(UseQ == kTRUE) fval = fval / jetsS->At(0)->Pt();
106	>	fval = fval / Q;
107
108		return fval;
109		}
110
111	<	Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, bool UseQ, Double_t Y){
111	>	Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, double Q, double Y){
112		if(calos->GetEntries() <= 0) return 0.0;
113
114		Double_t fval = 0.0;
#	Line 93 | Line 120 | Double_t JetTools::NJettiness(const Calo
120		for(int j=0;j<int(jets->GetEntries());j++){
121		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
122		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-calos->At(i)->Eta()))
123	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi()))));
123	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi())))));
124		}
125		fval = fval + fvalpart;
126		}
127
128	<	if(UseQ == kTRUE) fval = fval / calos->At(0)->Pt();
128	>	fval = fval / Q;
129
130		return fval;
131		}
#	Line 184 | Line 211 | Double_t JetTools::CosineOmega(const Par
211
212		//Transverse Higgs mass
213		Double_t JetTools::MtHiggs(const ParticleOArr * leptons,
214	<	const Met *met, int nsel){
214	>	const Met *met, double metFraction[2], int nsel){
215		if(leptons->Entries() < 2) return -999.0;
216
217		double mtHiggs = -999.0;
#	Line 230 | Line 257 | Double_t JetTools::MtHiggs(const Particl
257		mll   = dilepton->Mass();
258		mnu   = 0.0;
259		}
260	<	else if(nsel == 4){ // Use the formula from hep-ph:1006.4998
260	>	else if(nsel == 4){ // Use of Mt mass and replacing mnu using the met optimal
261	>	enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
262	>	enenn = TMath::Sqrt(met->Pt() *met->Pt()  + 0.0*0.0);
263	>	enex  = dilepton->Px() + met->Px();
264	>	eney  = dilepton->Py() + met->Py();
265	>	mll   = dilepton->Mass();
266	>	double metAuxPx[2] = {met->Px() * metFraction[0],
267	>	met->Px() * (1.0 - metFraction[0])};
268	>	double metAuxPy[2] = {met->Py() * metFraction[1],
269	>	met->Py() * (1.0 - metFraction[1])};
270	>	double ene = TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
271	>	TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
272	>	double px = metAuxPx[0] + metAuxPx[1];
273	>	double py = metAuxPy[0] + metAuxPy[1];
274	>	mnu = TMath::Sqrt(ene*ene - px*px - py*py);
275	>	}
276	>	else if(nsel == 5){ // Using the optimal met value
277	>	double metAuxPx[2] = {met->Px() * metFraction[0],
278	>	met->Px() * (1.0 - metFraction[0])};
279	>	double metAuxPy[2] = {met->Py() * metFraction[1],
280	>	met->Py() * (1.0 - metFraction[1])};
281	>	double ene = leptons->At(0)->Pt() + leptons->At(1)->Pt() +
282	>	TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
283	>	TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
284	>	double px = leptons->At(0)->Px() + leptons->At(1)->Px() +
285	>	metAuxPx[0] + metAuxPx[1];
286	>	double py = leptons->At(0)->Py() + leptons->At(1)->Py() +
287	>	metAuxPy[0] + metAuxPy[1];
288	>	mtHiggs = ene*ene - px*px - py*py;
289	>	}
290	>	else if(nsel == 6){ // Use the formula from hep-ph:1006.4998
291		mtHiggs = 2*leptons->At(0)->Pt()*leptons->At(0)->Pt() + 2*leptons->At(1)->Pt()*leptons->At(1)->Pt() + 3 * (
292		leptons->At(0)->Pt()*leptons->At(1)->Pt() + met->Pt()*(leptons->At(0)->Pt()+leptons->At(1)->Pt())
293		- met->Px()*dilepton->Px() - met->Py()*dilepton->Py()
294		- leptons->At(0)->Px()*leptons->At(1)->Px() - leptons->At(0)->Py()*leptons->At(1)->Py());
295		}
296	+	else if(nsel == 7){ // Use of M mass and mnu == 0
297	+	double deltaPhiDileptonMet = fabs(MathUtils::DeltaPhi(dilepton->Phi(),
298	+	met->Phi()));
299	+	mtHiggs = 2.0*dilepton->Pt()*met->Pt()*(1.0 - cos(deltaPhiDileptonMet));
300	+	}
301
302	<	if(nsel >= 0 && nsel <= 3){
302	>	if(nsel >= 0 && nsel <= 4){
303		mtHiggs = mll*mll + mnu*mnu + 2.0*(enell*enenn - enex*enex - eney*eney);
304		}
305	+
306		if(mtHiggs <= 0) mtHiggs = 0.0;
307		else             mtHiggs = TMath::Sqrt(mtHiggs);
308
#	Line 247 | Line 310 | Double_t JetTools::MtHiggs(const Particl
310
311		return mtHiggs;
312		}
313	+
314	+	Double_t JetTools::Beta(const TrackCol *tracks, Jet *jet, const Vertex *vertex, Double_t  delta_z, Double_t delta_cone){
315	+
316	+	if(tracks->GetEntries() <= 0) return 1.0;
317	+
318	+	double Pt_jets_X = 0. ;
319	+	double Pt_jets_Y = 0. ;
320	+	double Pt_jets_X_tot = 0. ;
321	+	double Pt_jets_Y_tot = 0. ;
322	+
323	+	for(int i=0;i<int(tracks->GetEntries());i++){
324	+	if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
325	+	Pt_jets_X_tot += tracks->At(i)->Px();
326	+	Pt_jets_Y_tot += tracks->At(i)->Py();
327	+	double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertex));
328	+	if(pDz < delta_z){
329	+	Pt_jets_X += tracks->At(i)->Px();
330	+	Pt_jets_Y += tracks->At(i)->Py();
331	+	}
332	+	}
333	+	}
334	+
335	+	if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
336	+	return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);
337	+
338	+	return 1.0;
339	+	}
340	+
341	+
342	+	Double_t JetTools::Beta(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){
343	+	double Pt_jets= 0. ;
344	+	double Pt_jetsTot = 0. ;
345	+
346	+	for(UInt_t i=0;i<jet->NPFCands();i++){
347	+	if(jet->PFCand(i)->TrackerTrk()){
348	+	Pt_jetsTot += jet->PFCand(i)->TrackerTrk()->Pt();
349	+	double pDz = TMath::Abs(jet->PFCand(i)->TrackerTrk()->DzCorrected(*vertex));
350	+	if(pDz < delta_z){
351	+	Pt_jets += jet->PFCand(i)->TrackerTrk()->Pt();
352	+	}
353	+	}
354	+	}
355	+
356	+	Double_t beta = 0.;
357	+	if (Pt_jetsTot > 0)
358	+	beta = Pt_jets/Pt_jetsTot;
359	+
360	+	return beta;
361	+	}
362	+
363	+	Double_t JetTools::Beta2(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){
364	+	double Pt_jets= 0. ;
365	+	double Pt_jetsTot = 0. ;
366	+
367	+	for(UInt_t i=0;i<jet->NPFCands();i++){
368	+	if(jet->PFCand(i)->BestTrk()){
369	+	Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
370	+	double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
371	+	if(pDz < delta_z){
372	+	Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
373	+	}
374	+	}
375	+	}
376	+
377	+	Double_t beta = 1.0;
378	+	if (Pt_jetsTot > 0)
379	+	beta = Pt_jets/Pt_jetsTot;
380	+	return beta;
381	+	}
382	+
383	+
384	+	Bool_t  JetTools::PassBetaVertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
385	+
386	+	Bool_t passBetaCut = kTRUE;
387	+	if(vertices->GetEntries() > 0) {
388	+	Double_t Beta = JetTools::Beta(jet, referenceVertex, 0.2);
389	+	Double_t Beta_other = 0.0;
390	+	for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
391	+	if (referenceVertex == vertices->At(nv)) continue;
392	+	Double_t BetaAux = JetTools::Beta(jet, vertices->At(nv), 0.2);
393	+	if(BetaAux > Beta_other) Beta_other = BetaAux;
394	+	}
395	+	if(Beta_other > Beta) passBetaCut = kFALSE;
396	+	}
397	+
398	+	return passBetaCut;
399	+
400	+	}
401	+
402	+	Bool_t  JetTools::PassBeta2VertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
403	+
404	+	Bool_t passBetaCut = kTRUE;
405	+	if(vertices->GetEntries() > 0) {
406	+	Double_t Beta = JetTools::Beta2(jet, referenceVertex, 0.2);
407	+	Double_t Beta_other = 0.0;
408	+	for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
409	+	if (referenceVertex == vertices->At(nv)) continue;
410	+	Double_t BetaAux = JetTools::Beta2(jet, vertices->At(nv), 0.2);
411	+	if(BetaAux > Beta_other) Beta_other = BetaAux;
412	+	}
413	+	if(Beta_other > Beta) passBetaCut = kFALSE;
414	+	}
415	+
416	+	return passBetaCut;
417	+
418	+	}
419	+
420	+
421	+	Int_t JetTools::MaxBetaVertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
422	+
423	+	Int_t vertexIndex = -1;
424	+	double beta = -0.1;
425	+	for (UInt_t v=0; v < vertices->GetEntries(); v++){
426	+	Double_t betaTmp = JetTools::Beta(jet, vertices->At(v), delta_z);
427	+	if (betaTmp > beta) {
428	+	beta = betaTmp;
429	+	vertexIndex = v;
430	+	}
431	+	}
432	+	return vertexIndex;
433	+
434	+	}
435	+
436	+	Int_t JetTools::MaxBeta2VertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
437	+
438	+	Int_t vertexIndex = -1;
439	+	double beta = -0.1;
440	+	for (UInt_t v=0; v < vertices->GetEntries(); v++){
441	+	Double_t betaTmp = JetTools::Beta2(jet, vertices->At(v), delta_z);
442	+	if (betaTmp > beta) {
443	+	beta = betaTmp;
444	+	vertexIndex = v;
445	+	}
446	+	}
447	+	return vertexIndex;
448	+
449	+	}
450	+
451	+
452	+	Int_t JetTools::JetToPVAssociation(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
453	+
454	+	std::vector<float> verticesPt2(vertices->GetEntries());
455	+	for(UInt_t i=0;i<jet->NPFCands();i++){
456	+	if(jet->PFCand(i)->BestTrk()){
457	+	double minDZ = delta_z;
458	+	int trackVertexIndex = -1;
459	+	for (UInt_t v=0; v < vertices->GetEntries(); v++){
460	+	if (minDZ > TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)))) {
461	+	minDZ = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)));
462	+	trackVertexIndex = v;
463	+	}
464	+	}
465	+	if (trackVertexIndex < 0) continue;
466	+	verticesPt2[trackVertexIndex]+= jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
467	+	}
468	+	}
469	+
470	+	Int_t vertexIndex = 0;
471	+	float pt2Max = 0;
472	+	for (uint i=0; i < verticesPt2.size(); ++i){
473	+	if (pt2Max < verticesPt2[i]) {
474	+	pt2Max = verticesPt2[i];
475	+	vertexIndex = i;
476	+	}
477	+	}
478	+	return vertexIndex;
479	+	}
480	+	const PFCandidate* JetTools::leadCand(const PFJet *iJet,int iPFType,bool i2nd) {
481	+	int lCount = 0;
482	+	const PFCandidate *lCand = 0;
483	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
484	+	lCand = iJet->PFCand(i0);
485	+	if(iPFType != -1 && lCand->PFType() != iPFType) continue;
486	+	if(lCount == 0 && !i2nd) break;
487	+	if(lCount >  0)          break;
488	+	lCount++;
489	+	}
490	+	return lCand;
491	+	}
492	+	Double_t JetTools::impactParameter(const PFJet *iJet,const Vertex *iVertex,bool iDZ) {
493	+	double lDZCorr = -1000;
494	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
495	+	const PFCandidate *pCand = iJet->PFCand(i0);
496	+	if(pCand->TrackerTrk() == 0) continue;
497	+	//if(pCand->Pt() < 1.) continue; => previous iterations
498	+	if(iDZ)  lDZCorr = fabs(pCand->TrackerTrk()->DzCorrected(*iVertex));
499	+	if(!iDZ) lDZCorr = fabs(pCand->TrackerTrk()->D0Corrected(*iVertex));
500	+	break;
501	+	}
502	+	return lDZCorr;
503	+	}
504	+	Double_t JetTools::dRMean(const PFJet *iJet,int iPFType) {
505	+	double lDRMean = 0;
506	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
507	+	const PFCandidate *pCand = iJet->PFCand(i0);
508	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
509	+	double pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
510	+	lDRMean    += pDR*(pCand->Pt())/iJet->RawMom().Pt();
511	+	}
512	+	return lDRMean;
513	+	}
514	+	Double_t JetTools::dR2Mean(const PFJet *iJet,int iPFType) {
515	+	double lDR2Mean = 0;
516	+	double lSumPt2 = 0;
517	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
518	+	const PFCandidate *pCand = iJet->PFCand(i0);
519	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
520	+	lSumPt2   += pCand->Pt()*pCand->Pt();
521	+	double pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
522	+	lDR2Mean    += pDR*pDR*(pCand->Pt()*pCand->Pt());
523	+	}
524	+	lDR2Mean/=lSumPt2;
525	+	return lDR2Mean;
526	+	}
527	+	Double_t JetTools::frac(const PFJet *iJet,Double_t iDRMax,Double_t iDRMin,Int_t iPFType) {
528	+	double lFrac = 0;
529	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
530	+	const PFCandidate *pCand = iJet->PFCand(i0);
531	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
532	+	Double_t pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
533	+	if(pDR > iDRMax) continue;
534	+	if(pDR < iDRMax-0.1) continue;
535	+	lFrac += pCand->Pt()/iJet->RawMom().Pt();
536	+	}
537	+	return lFrac;
538	+	}
539	+	Double_t JetTools::betaStar(const PFJet *iJet,const Vertex *iVertex,const VertexCol* iVertices,Double_t iDZCut) {
540	+	Double_t lTotal = 0;
541	+	Double_t lPileup = 0;
542	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
543	+	const PFCandidate* pPF   = iJet->PFCand(i0);
544	+	const Track* pTrack      = pPF->TrackerTrk();
545	+	//if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW
546	+	if(pTrack == 0) continue;
547	+	lTotal += pTrack->Pt();
548	+	double pDZPV  = fabs(pTrack->DzCorrected(*iVertex));
549	+	double pDZMin = pDZPV;
550	+	for(unsigned int i1 = 0; i1 < iVertices->GetEntries(); i1++) {
551	+	const Vertex *pV = iVertices->At(i1);
552	+	if(pV->Ndof() < 4 ||
553	+	(pV->Position() - iVertex->Position()).R() < 0.02 ) continue;
554	+	pDZMin = TMath::Min(pDZMin,fabs(pTrack->DzCorrected(*pV)));
555	+	}
556	+	if(pDZPV > 0.2 && pDZMin < 0.2) lPileup += pTrack->Pt();
557	+	}
558	+	if(lTotal == 0) lTotal = 1;
559	+	return lPileup/(lTotal);
560	+	}
561	+	Double_t JetTools::betaStarClassic(const PFJet *iJet,const Vertex *iVertex,const VertexCol* iVertices) {
562	+	Double_t lTotal = 0;
563	+	Double_t lPileup = 0;
564	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
565	+	const PFCandidate* pPF   = iJet->PFCand(i0);
566	+	const Track* pTrack      = pPF->TrackerTrk();
567	+	//if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW
568	+	if(pTrack == 0) continue;
569	+	lTotal += pTrack->Pt();
570	+	bool isPV    = iVertex->HasTrack(pPF->TrackerTrk());
571	+	bool isOtherV = false;
572	+	for(unsigned int i1 = 0; i1 < iVertices->GetEntries(); i1++) {
573	+	const Vertex *pV = iVertices->At(i1);
574	+	if(isOtherV || isPV) continue;
575	+	if(pV->Ndof() < 4 ||
576	+	(pV->Position() - iVertex->Position()).R() < 0.02 ) continue;
577	+	isOtherV    = pV->HasTrack(pPF->TrackerTrk());
578	+	}
579	+	if(!isPV && isOtherV) lPileup += pTrack->Pt();
580	+	}
581	+	if(lTotal == 0) lTotal = 1;
582	+	return lPileup/(lTotal);
583	+	}
584	+	Bool_t  JetTools::passPFLooseId(const PFJet *iJet) {
585	+	if(iJet->RawMom().E()                              == 0)       return false;
586	+	if(iJet->NeutralHadronEnergy()/iJet->RawMom().E()  >  0.99)    return false;
587	+	if(iJet->NeutralEmEnergy()/iJet->RawMom().E()      >  0.99)    return false;
588	+	if(iJet->NConstituents()                           <  2)       return false;
589	+	if(iJet->ChargedHadronEnergy()/iJet->RawMom().E()  <= 0     && fabs(iJet->Eta()) < 2.4 ) return false;
590	+	if(iJet->ChargedEmEnergy()/iJet->RawMom().E()      >  0.99  && fabs(iJet->Eta()) < 2.4 ) return false;
591	+	if(iJet->ChargedMultiplicity()                     < 1      && fabs(iJet->Eta()) < 2.4 ) return false;
592	+	//if(fabs(iJet->Eta())                               > 4.99) return false;
593	+	return true;
594	+	}
595	+	//Jet Width Variables
596	+	double JetTools::W(const PFJet *iJet,int iPFType,int iType) {
597	+	double lPtD    = 0;
598	+	double lSumPt  = 0;
599	+	double lSumPt2 = 0;
600	+	TMatrixDSym lCovMatrix(2); lCovMatrix = 0.;
601	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
602	+	const PFCandidate *pCand = iJet->PFCand(i0);
603	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
604	+	double pDEta = iJet->Eta() - pCand->Eta();
605	+	double pDPhi = fabs(iJet->Phi()-pCand->Phi()); if(pDPhi > 2.*TMath::Pi() - pDPhi) pDPhi =  2.*TMath::Pi() - pDPhi;
606	+	lCovMatrix(0,0) += pCand->Pt()*pCand->Pt()*pDEta*pDEta;
607	+	lCovMatrix(0,1) += pCand->Pt()*pCand->Pt()*pDEta*pDPhi;
608	+	lCovMatrix(1,1) += pCand->Pt()*pCand->Pt()*pDPhi*pDPhi;
609	+	lPtD            += pCand->Pt()*pCand->Pt();
610	+	lSumPt          += pCand->Pt();
611	+	lSumPt2         += pCand->Pt()*pCand->Pt();
612	+	}
613	+	lCovMatrix(0,0) /= lSumPt2;
614	+	lCovMatrix(0,1) /= lSumPt2;
615	+	lCovMatrix(1,1) /= lSumPt2;
616	+	lCovMatrix(1,0)  = lCovMatrix(0,1);
617	+	lPtD             = sqrt(lPtD);
618	+	lPtD            /= lSumPt;
619	+	double lEtaW     = sqrt(lCovMatrix(0,0));
620	+	double lPhiW     = sqrt(lCovMatrix(1,1));
621	+	double lJetW     = 0.5*(lEtaW+lPhiW);
622	+	TVectorD lEigVals(2); lEigVals = TMatrixDSymEigen(lCovMatrix).GetEigenValues();
623	+	double lMajW     = sqrt(fabs(lEigVals(0)));
624	+	double lMinW     = sqrt(fabs(lEigVals(1)));
625	+
626	+	if(iType == 1) return lMajW;
627	+	if(iType == 2) return lMinW;
628	+	if(iType == 3) return lEtaW;
629	+	if(iType == 4) return lPhiW;
630	+	if(iType == 5) return lJetW;
631	+	return lPtD; //ptRMS
632	+	}
633	+	/*
634	+	double JetTools::genFrac(const PFJet *iJet) {
635	+	double lTrueFrac = 0;
636	+	for(UInt_t i0 = 0; i0 < fParticles->GetEntries(); i0++) {
637	+	const MCParticle *p = fParticles->At(i0);
638	+	if(p->Status() != 1) continue;
639	+	double pDEta = iJet->Eta() - p->Eta();
640	+	double pDPhi = fabs(iJet->Phi()-p->Phi()); if(pDPhi > 2.*TMath::Pi() - pDPhi) pDPhi =  2.*TMath::Pi() - pDPhi;
641	+	double pDR   = sqrt(pDEta*pDEta + pDPhi*pDPhi);
642	+	if(pDR > 0.5) continue;
643	+	lTrueFrac += p->Pt();
644	+	}
645	+	lTrueFrac/=iJet->Pt();
646	+	return lTrueFrac;
647	+	}
648	+	*/
649	+	/*
650	+	double* JetTools::subStructure(const PFJet *iJet) {
651	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
652	+	const PFCandidate *pCand = iJet->PFCand(i0);
653	+	//Fast Jet
654	+	FJparticles.push_back( fastjet::PseudoJet( pCand->Px(),pCand->Py(),pCand->Pz(),pCand->Energy() ) );
655	+	}
656	+	//ReCluster
657	+	fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, 0.5);
658	+	int activeAreaRepeats = 1;
659	+	double ghostArea = 0.01;
660	+	double ghostEtaMax = 5.0;
661	+	fastjet::ActiveAreaSpec fjActiveArea(ghostEtaMax,activeAreaRepeats,ghostArea);
662	+	fjActiveArea.set_fj2_placement(true);
663	+
664	+	fastjet::AreaDefinition fjAreaDefinition(fastjet::active_area_explicit_ghosts, fjActiveArea );
665	+	fastjet::ClusterSequenceArea thisClustering(FJparticles, jetDef, fjAreaDefinition);
666	+	std::vector<fastjet::PseudoJet> out_jets = sorted_by_pt(thisClustering.inclusive_jets(0.0));
667	+	//std::cout << "===>  Size " << FJparticles.size() << " -- " << out_jets.size() << std::endl;
668	+	// Adding substructure
669	+	// define groomers
670	+	fastjet::Filter trimmer( fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.03)));
671	+	fastjet::Filter filter( fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, 0.3), fastjet::SelectorNHardest(3)));
672	+	fastjet::Pruner pruner(fastjet::cambridge_algorithm, 0.1, 0.5);
673	+
674	+	std::vector<fastjet::Transformer const *> transformers;
675	+	transformers.push_back(&trimmer);
676	+	transformers.push_back(&filter);
677	+	transformers.push_back(&pruner);
678	+
679	+	// define n-subjettiness
680	+	NsubParameters paraNsub = NsubParameters(1.0, 0.8);
681	+	Nsubjettiness routine(nsub_kt_axes, paraNsub);
682	+	for (unsigned i0 = 0; i0 < out_jets.size(); i0++) {
683	+	int i1 = -1;
684	+	if(out_jets.at(i0).pt()  < 1.) continue;
685	+	//std::cout << " ---> " << out_jets.at(i0).pt() << std::endl;
686	+	for ( std::vector<fastjet::Transformer const *>::const_iterator itransf = transformers.begin(), itransfEnd = transformers.end(); itransf != itransfEnd; ++itransf ) {
687	+	i1++;
688	+	fastjet::PseudoJet transformedJet = out_jets.at(i0);
689	+	transformedJet = (**itransf)(transformedJet);
690	+	TLorentzVector jetcorr(transformedJet.px() * jec,transformedJet.py() * jec,transformedJet.pz() * jec,transformedJet.e() * jec);
691	+	if(i1 == 0) internalId_.trimmass_   = jetcorr.M();
692	+	if(i1 == 0) internalId_.trimarea_   = transformedJet.area();
693	+	if(i1 == 1) internalId_.filtermass_ = jetcorr.M();
694	+	if(i1 == 1) internalId_.filterarea_ = transformedJet.area();
695	+	if(i1 == 2) internalId_.prunedmass_ = jetcorr.M();
696	+	if(i1 == 2) internalId_.prunedarea_ = transformedJet.area();
697	+	if (transformedJet.constituents().size() > 1 && i1 == 2 ) {
698	+	std::vector<fastjet::PseudoJet> subjets = transformedJet.associated_cluster_sequence()->exclusive_subjets(transformedJet,2);
699	+	internalId_.nsubjets_     = subjets.size();
700	+	internalId_.massdrop_     = subjets.at(0).m()/transformedJet.m();
701	+	internalId_.massdropcorr_ = subjets.at(0).m()/internalId_.prunedmass_;
702	+	}
703	+	}
704	+	internalId_.tau1_ = routine.getTau(1, out_jets.at(i0).constituents());
705	+	internalId_.tau2_ = routine.getTau(2, out_jets.at(i0).constituents());
706	+	internalId_.tau3_ = routine.getTau(3, out_jets.at(i0).constituents());
707	+	internalId_.tau4_ = routine.getTau(4, out_jets.at(i0).constituents());
708	+	internalId_.tau2tau1_ = internalId_.tau2_/internalId_.tau1_;
709	+	}
710	+	}
711	+	}
712	+	*/

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines