ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetTools.cc
Revision: 1.29
Committed: Mon May 14 09:42:27 2012 UTC (12 years, 11 months ago) by pharris
Content type: text/plain
Branch: MAIN
Changes since 1.28: +13 -0 lines
Log Message: 
Added final MVA

File Contents

# User Rev Content
1 ceballos 1.1 #include "MitPhysics/Utils/interface/JetTools.h"
2 mzanetti 1.13 #include <algorithm>
3     #include <vector>
4 ceballos 1.1
5     ClassImp(mithep::JetTools)
6    
7     using namespace mithep;
8 mzanetti 1.13
9 ceballos 1.1
10     JetTools::JetTools()
11     {
12     // Constructor
13     }
14    
15     JetTools::~JetTools()
16     {
17     // Destructor.
18     }
19    
20     //Remember to remove the signal from particles before inputting into the function
21 mzanetti 1.6 Double_t JetTools::NJettiness(const ParticleOArr *particles, const JetOArr *jets, double Q, double Y){
22 ceballos 1.2 if(particles->GetEntries() <= 0) return 0.0;
23    
24 ceballos 1.1 Double_t fval = 0.0;
25     Double_t fvalpart;
26    
27     for(int i=0;i<int(particles->GetEntries());i++){
28     fvalpart = (particles->At(i)->Pt()) * TMath::Exp(-TMath::Abs(particles->At(i)->Eta()-Y));
29    
30     for(int j=0;j<int(jets->GetEntries());j++){
31     fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
32     (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-particles->At(i)->Eta()))
33 sixie 1.15 - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi())))));
34 ceballos 1.1 }
35     fval = fval + fvalpart;
36     }
37 ceballos 1.2
38 mzanetti 1.6 fval = fval / Q;
39 ceballos 1.2
40 ceballos 1.1 return fval;
41     }
42    
43 mzanetti 1.6 Double_t JetTools::NJettiness(const PFCandidateOArr *pfCandidates, const JetOArr *jets, double Q, double Y){
44     if(pfCandidates->GetEntries() <= 0) return 0.0;
45    
46     Double_t fval = 0.0;
47     Double_t fvalpart;
48    
49     for(int i=0;i<int(pfCandidates->GetEntries());i++){
50     fvalpart = (pfCandidates->At(i)->Pt()) * TMath::Exp(-TMath::Abs(pfCandidates->At(i)->Eta()-Y));
51    
52     for(int j=0;j<int(jets->GetEntries());j++){
53     fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
54     (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-pfCandidates->At(i)->Eta()))
55 sixie 1.15 - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi())))));
56 mzanetti 1.6 }
57     fval = fval + fvalpart;
58     }
59    
60     fval = fval / Q;
61    
62     return fval;
63     }
64    
65     Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, double Q, double Y){
66 ceballos 1.2 if(tracks->GetEntries() <= 0) return 0.0;
67    
68 ceballos 1.1 Double_t fval = 0.0;
69     Double_t fvalpart;
70    
71     for(int i=0;i<int(tracks->GetEntries());i++){
72     fvalpart = (tracks->At(i)->Pt()) * TMath::Exp(-TMath::Abs(tracks->At(i)->Eta()-Y));
73    
74     for(int j=0;j<int(jets->GetEntries());j++){
75     fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
76     (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-tracks->At(i)->Eta()))
77 sixie 1.15 - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi())))));
78 ceballos 1.1 }
79     fval = fval + fvalpart;
80     }
81 ceballos 1.2
82 mzanetti 1.6 fval = fval / Q;
83 ceballos 1.2
84 ceballos 1.1 return fval;
85     }
86    
87 mzanetti 1.6 Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, double Q, double Y){
88 ceballos 1.2 if(jetsS->GetEntries() <= 0) return 0.0;
89    
90 ceballos 1.1 Double_t fval = 0.0;
91     Double_t fvalpart;
92    
93     for(int i=0;i<int(jetsS->GetEntries());i++){
94     fvalpart = (jetsS->At(i)->Pt()) * TMath::Exp(-TMath::Abs(jetsS->At(i)->Eta()-Y));
95    
96     for(int j=0;j<int(jets->GetEntries());j++){
97     fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
98     (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-jetsS->At(i)->Eta()))
99 sixie 1.15 - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi())))));
100 ceballos 1.1 }
101     fval = fval + fvalpart;
102     }
103 ceballos 1.2
104 mzanetti 1.6 fval = fval / Q;
105 ceballos 1.2
106     return fval;
107     }
108    
109 mzanetti 1.6 Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, double Q, double Y){
110 ceballos 1.2 if(calos->GetEntries() <= 0) return 0.0;
111    
112     Double_t fval = 0.0;
113     Double_t fvalpart;
114    
115     for(int i=0;i<int(calos->GetEntries());i++){
116     fvalpart = (calos->At(i)->Pt()) * TMath::Exp(-TMath::Abs(calos->At(i)->Eta()-Y));
117    
118     for(int j=0;j<int(jets->GetEntries());j++){
119     fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
120     (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-calos->At(i)->Eta()))
121 sixie 1.15 - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi())))));
122 ceballos 1.2 }
123     fval = fval + fvalpart;
124     }
125    
126 mzanetti 1.6 fval = fval / Q;
127 ceballos 1.2
128 ceballos 1.1 return fval;
129     }
130    
131     //M_r
132     Double_t JetTools::M_r(const ParticleOArr *particles){
133    
134     if(particles->GetEntries() < 2) return -999.;
135    
136     Double_t E0 = particles->At(0)->E();
137     Double_t E1 = particles->At(1)->E();
138     Double_t Pz0 = particles->At(0)->Pz();
139     Double_t Pz1 = particles->At(1)->Pz();
140    
141     Double_t den = TMath::Power(Pz0-Pz1, 2) - TMath::Power(E0-E1,2);
142     if(den <= 0) return -100.;
143    
144     return 2.0*TMath::Sqrt(TMath::Power(E0*Pz1 - E1*Pz0, 2)/den);
145     }
146    
147     //Beta_r
148     Double_t JetTools::Beta_r(const ParticleOArr *particles){
149    
150     if(particles->GetEntries() < 2) return -999.;
151    
152     Double_t E0 = particles->At(0)->E();
153     Double_t E1 = particles->At(1)->E();
154     Double_t Pz0 = particles->At(0)->Pz();
155     Double_t Pz1 = particles->At(1)->Pz();
156    
157     return (E0-E1)/(Pz0-Pz1);
158     }
159    
160     //M_r_t
161     Double_t JetTools::M_r_t(const ParticleOArr *particles, const Met *met){
162    
163     if(particles->GetEntries() < 2) return -999.;
164    
165     Double_t Pt0 = particles->At(0)->Pt();
166     Double_t Pt1 = particles->At(1)->Pt();
167     Double_t etmiss = met->Pt();
168    
169     Double_t Px0 = particles->At(0)->Px();
170     Double_t Px1 = particles->At(1)->Px();
171     Double_t metx = met->Px();
172     Double_t Py0 = particles->At(0)->Py();
173     Double_t Py1 = particles->At(1)->Py();
174     Double_t mety = met->Py();
175    
176     return TMath::Sqrt(0.5*etmiss*(Pt0 + Pt1) - 0.5*(metx*(Px0 + Px1) + mety*(Py0 + Py1)));
177     }
178    
179     //Razor
180     Double_t JetTools::Razor(const ParticleOArr *particles, const Met *met){
181     if(particles->GetEntries() < 2) return -999.;
182    
183     Double_t mr = M_r(particles);
184     Double_t mrt = M_r_t(particles,met);
185    
186     if(mr != 0) return mrt/mr;
187    
188     return -999.;
189     }
190    
191     //Cosine Omega
192 ceballos 1.3 Double_t JetTools::CosineOmega(const Particle *particles0, const Particle *particles1){
193 ceballos 1.1
194 ceballos 1.3 TLorentzVector v_L1(particles0->Px(),particles0->Py(),particles0->Pz(),particles0->E());
195     TLorentzVector v_L2(particles1->Px(),particles1->Py(),particles1->Pz(),particles1->E());
196 ceballos 1.1
197     Double_t beta = (v_L1.P()-v_L2.P())/(v_L1.Pz()-v_L2.Pz());
198    
199     TVector3 B;
200     B.SetXYZ(0.0,0.0,-1.0*beta);
201    
202     v_L1.Boost(B);
203     v_L2.Boost(B);
204    
205     Double_t cosomega = v_L1.Vect().Dot(v_L2.Vect())/(v_L1.P()*v_L2.P());
206    
207     return cosomega;
208     }
209    
210     //Transverse Higgs mass
211 ceballos 1.3 Double_t JetTools::MtHiggs(const ParticleOArr * leptons,
212 ceballos 1.4 const Met *met, double metFraction[2], int nsel){
213 ceballos 1.3 if(leptons->Entries() < 2) return -999.0;
214    
215 ceballos 1.1 double mtHiggs = -999.0;
216 ceballos 1.3 double enell = 0.0;
217     double enenn = 0.0;
218     double enex = 0.0;
219     double eney = 0.0;
220     double mll = 0.0;
221     double mnu = 0.0;
222     CompositeParticle *dilepton = new CompositeParticle();
223     dilepton->AddDaughter(leptons->At(0));
224     dilepton->AddDaughter(leptons->At(1));
225 ceballos 1.1
226     if (nsel == 0){ // Use of Mt mass and mnu == mll
227     enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
228     enenn = TMath::Sqrt(met->Pt() *met->Pt() + dilepton->Mt()*dilepton->Mt());
229     enex = dilepton->Px() + met->Px();
230     eney = dilepton->Py() + met->Py();
231     mll = dilepton->Mass();
232     mnu = mll;
233     }
234     else if(nsel == 1){ // Use of Mt mass and mnu == 0
235     enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
236     enenn = TMath::Sqrt(met->Pt() *met->Pt() + 0.0*0.0);
237     enex = dilepton->Px() + met->Px();
238     eney = dilepton->Py() + met->Py();
239     mll = dilepton->Mass();
240     mnu = 0.0;
241     }
242     else if(nsel == 2){ // Use of M mass and mnu == mll
243     enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mass()*dilepton->Mass());
244     enenn = TMath::Sqrt(met->Pt() *met->Pt() + dilepton->Mass()*dilepton->Mass());
245     enex = dilepton->Px() + met->Px();
246     eney = dilepton->Py() + met->Py();
247     mll = dilepton->Mass();
248     mnu = mll;
249     }
250     else if(nsel == 3){ // Use of M mass and mnu == 0
251     enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mass()*dilepton->Mass());
252     enenn = TMath::Sqrt(met->Pt() *met->Pt() + 0.0*0.0);
253     enex = dilepton->Px() + met->Px();
254     eney = dilepton->Py() + met->Py();
255     mll = dilepton->Mass();
256     mnu = 0.0;
257     }
258 ceballos 1.4 else if(nsel == 4){ // Use of Mt mass and replacing mnu using the met optimal
259     enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
260     enenn = TMath::Sqrt(met->Pt() *met->Pt() + 0.0*0.0);
261     enex = dilepton->Px() + met->Px();
262     eney = dilepton->Py() + met->Py();
263     mll = dilepton->Mass();
264     double metAuxPx[2] = {met->Px() * metFraction[0],
265     met->Px() * (1.0 - metFraction[0])};
266     double metAuxPy[2] = {met->Py() * metFraction[1],
267     met->Py() * (1.0 - metFraction[1])};
268     double ene = TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
269     TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
270     double px = metAuxPx[0] + metAuxPx[1];
271     double py = metAuxPy[0] + metAuxPy[1];
272     mnu = TMath::Sqrt(ene*ene - px*px - py*py);
273     }
274     else if(nsel == 5){ // Using the optimal met value
275     double metAuxPx[2] = {met->Px() * metFraction[0],
276     met->Px() * (1.0 - metFraction[0])};
277     double metAuxPy[2] = {met->Py() * metFraction[1],
278     met->Py() * (1.0 - metFraction[1])};
279     double ene = leptons->At(0)->Pt() + leptons->At(1)->Pt() +
280     TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
281     TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
282     double px = leptons->At(0)->Px() + leptons->At(1)->Px() +
283     metAuxPx[0] + metAuxPx[1];
284     double py = leptons->At(0)->Py() + leptons->At(1)->Py() +
285     metAuxPy[0] + metAuxPy[1];
286     mtHiggs = ene*ene - px*px - py*py;
287     }
288     else if(nsel == 6){ // Use the formula from hep-ph:1006.4998
289 ceballos 1.3 mtHiggs = 2*leptons->At(0)->Pt()*leptons->At(0)->Pt() + 2*leptons->At(1)->Pt()*leptons->At(1)->Pt() + 3 * (
290     leptons->At(0)->Pt()*leptons->At(1)->Pt() + met->Pt()*(leptons->At(0)->Pt()+leptons->At(1)->Pt())
291     - met->Px()*dilepton->Px() - met->Py()*dilepton->Py()
292     - leptons->At(0)->Px()*leptons->At(1)->Px() - leptons->At(0)->Py()*leptons->At(1)->Py());
293 ceballos 1.1 }
294 ceballos 1.14 else if(nsel == 7){ // Use of M mass and mnu == 0
295 sixie 1.15 double deltaPhiDileptonMet = fabs(MathUtils::DeltaPhi(dilepton->Phi(),
296     met->Phi()));
297 ceballos 1.14 mtHiggs = 2.0*dilepton->Pt()*met->Pt()*(1.0 - cos(deltaPhiDileptonMet));
298 ceballos 1.11 }
299 ceballos 1.1
300 ceballos 1.4 if(nsel >= 0 && nsel <= 4){
301 ceballos 1.3 mtHiggs = mll*mll + mnu*mnu + 2.0*(enell*enenn - enex*enex - eney*eney);
302     }
303 ceballos 1.11
304 ceballos 1.1 if(mtHiggs <= 0) mtHiggs = 0.0;
305     else mtHiggs = TMath::Sqrt(mtHiggs);
306    
307 ceballos 1.3 delete dilepton;
308    
309 ceballos 1.1 return mtHiggs;
310     }
311 ceballos 1.5
312 ceballos 1.9 Double_t JetTools::Beta(const TrackCol *tracks, Jet *jet, const Vertex *vertex, Double_t delta_z, Double_t delta_cone){
313    
314     if(tracks->GetEntries() <= 0) return 1.0;
315 ceballos 1.5
316     double Pt_jets_X = 0. ;
317     double Pt_jets_Y = 0. ;
318     double Pt_jets_X_tot = 0. ;
319     double Pt_jets_Y_tot = 0. ;
320    
321     for(int i=0;i<int(tracks->GetEntries());i++){
322     if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
323     Pt_jets_X_tot += tracks->At(i)->Px();
324     Pt_jets_Y_tot += tracks->At(i)->Py();
325 ceballos 1.9 double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertex));
326 ceballos 1.5 if(pDz < delta_z){
327     Pt_jets_X += tracks->At(i)->Px();
328     Pt_jets_Y += tracks->At(i)->Py();
329     }
330     }
331     }
332    
333 ceballos 1.8 if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
334 ceballos 1.9 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);
335    
336     return 1.0;
337 ceballos 1.5 }
338    
339 ceballos 1.7
340 ceballos 1.9 Double_t JetTools::Beta(const PFJet *jet, const Vertex *vertex, Double_t delta_z){
341 mzanetti 1.13 double Pt_jets= 0. ;
342     double Pt_jetsTot = 0. ;
343 pharris 1.19
344 ceballos 1.7 for(UInt_t i=0;i<jet->NPFCands();i++){
345 pharris 1.21 if(jet->PFCand(i)->TrackerTrk()){
346     Pt_jetsTot += jet->PFCand(i)->TrackerTrk()->Pt();
347     double pDz = TMath::Abs(jet->PFCand(i)->TrackerTrk()->DzCorrected(*vertex));
348 ceballos 1.7 if(pDz < delta_z){
349 pharris 1.21 Pt_jets += jet->PFCand(i)->TrackerTrk()->Pt();
350 ceballos 1.7 }
351     }
352     }
353    
354 pharris 1.19 Double_t beta = 0.;
355 mzanetti 1.13 if (Pt_jetsTot > 0)
356     beta = Pt_jets/Pt_jetsTot;
357    
358     return beta;
359     }
360    
361     Double_t JetTools::Beta2(const PFJet *jet, const Vertex *vertex, Double_t delta_z){
362     double Pt_jets= 0. ;
363     double Pt_jetsTot = 0. ;
364    
365     for(UInt_t i=0;i<jet->NPFCands();i++){
366     if(jet->PFCand(i)->BestTrk()){
367     Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
368     double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
369     if(pDz < delta_z){
370     Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
371     }
372     }
373     }
374 ceballos 1.9
375 mzanetti 1.13 Double_t beta = 1.0;
376     if (Pt_jetsTot > 0)
377     beta = Pt_jets/Pt_jetsTot;
378     return beta;
379 ceballos 1.7 }
380 sixie 1.10
381    
382     Bool_t JetTools::PassBetaVertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
383    
384     Bool_t passBetaCut = kTRUE;
385     if(vertices->GetEntries() > 0) {
386     Double_t Beta = JetTools::Beta(jet, referenceVertex, 0.2);
387     Double_t Beta_other = 0.0;
388     for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
389     if (referenceVertex == vertices->At(nv)) continue;
390     Double_t BetaAux = JetTools::Beta(jet, vertices->At(nv), 0.2);
391     if(BetaAux > Beta_other) Beta_other = BetaAux;
392     }
393     if(Beta_other > Beta) passBetaCut = kFALSE;
394     }
395    
396     return passBetaCut;
397    
398     }
399 mzanetti 1.13
400     Bool_t JetTools::PassBeta2VertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
401    
402     Bool_t passBetaCut = kTRUE;
403     if(vertices->GetEntries() > 0) {
404     Double_t Beta = JetTools::Beta2(jet, referenceVertex, 0.2);
405     Double_t Beta_other = 0.0;
406     for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
407     if (referenceVertex == vertices->At(nv)) continue;
408     Double_t BetaAux = JetTools::Beta2(jet, vertices->At(nv), 0.2);
409     if(BetaAux > Beta_other) Beta_other = BetaAux;
410     }
411     if(Beta_other > Beta) passBetaCut = kFALSE;
412     }
413    
414     return passBetaCut;
415    
416     }
417    
418    
419     Int_t JetTools::MaxBetaVertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t delta_z=0.2){
420    
421     Int_t vertexIndex = -1;
422     double beta = -0.1;
423     for (UInt_t v=0; v < vertices->GetEntries(); v++){
424     Double_t betaTmp = JetTools::Beta(jet, vertices->At(v), delta_z);
425     if (betaTmp > beta) {
426     beta = betaTmp;
427     vertexIndex = v;
428     }
429     }
430     return vertexIndex;
431    
432     }
433    
434     Int_t JetTools::MaxBeta2VertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t delta_z=0.2){
435    
436     Int_t vertexIndex = -1;
437     double beta = -0.1;
438     for (UInt_t v=0; v < vertices->GetEntries(); v++){
439     Double_t betaTmp = JetTools::Beta2(jet, vertices->At(v), delta_z);
440     if (betaTmp > beta) {
441     beta = betaTmp;
442     vertexIndex = v;
443     }
444     }
445     return vertexIndex;
446    
447     }
448    
449    
450     Int_t JetTools::JetToPVAssociation(const PFJet *jet, const VertexCol *vertices, Double_t delta_z=0.2){
451    
452     std::vector<float> verticesPt2(vertices->GetEntries());
453     for(UInt_t i=0;i<jet->NPFCands();i++){
454     if(jet->PFCand(i)->BestTrk()){
455     double minDZ = delta_z;
456     int trackVertexIndex = -1;
457     for (UInt_t v=0; v < vertices->GetEntries(); v++){
458     if (minDZ > TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)))) {
459     minDZ = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)));
460     trackVertexIndex = v;
461     }
462     }
463     if (trackVertexIndex < 0) continue;
464     verticesPt2[trackVertexIndex]+= jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
465     }
466     }
467    
468     Int_t vertexIndex = 0;
469     float pt2Max = 0;
470     for (uint i=0; i < verticesPt2.size(); ++i){
471     if (pt2Max < verticesPt2[i]) {
472     pt2Max = verticesPt2[i];
473     vertexIndex = i;
474     }
475     }
476     return vertexIndex;
477     }
478 pharris 1.16 const PFCandidate* JetTools::leadCand(const PFJet *iJet,int iPFType,bool i2nd) {
479     int lCount = 0;
480     const PFCandidate *lCand = 0;
481     for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
482     lCand = iJet->PFCand(i0);
483     if(iPFType != -1 && lCand->PFType() != iPFType) continue;
484     if(lCount == 0 && !i2nd) break;
485     if(lCount > 0) break;
486     lCount++;
487     }
488     return lCand;
489     }
490     Double_t JetTools::impactParameter(const PFJet *iJet,const Vertex *iVertex,bool iDZ) {
491     double lDZCorr = -1000;
492     for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
493     const PFCandidate *pCand = iJet->PFCand(i0);
494 ceballos 1.22 if(pCand->TrackerTrk() == 0) continue;
495 pharris 1.19 //if(pCand->Pt() < 1.) continue; => previous iterations
496 pharris 1.21 if(iDZ) lDZCorr = fabs(pCand->TrackerTrk()->DzCorrected(*iVertex));
497     if(!iDZ) lDZCorr = fabs(pCand->TrackerTrk()->D0Corrected(*iVertex));
498 pharris 1.16 break;
499     }
500     return lDZCorr;
501     }
502     Double_t JetTools::dRMean(const PFJet *iJet,int iPFType) {
503     double lDRMean = 0;
504     for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
505     const PFCandidate *pCand = iJet->PFCand(i0);
506     if(iPFType != -1 && pCand->PFType() != iPFType) continue;
507     double pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
508 pharris 1.17 lDRMean += pDR*(pCand->Pt())/iJet->RawMom().Pt();
509 pharris 1.16 }
510     return lDRMean;
511     }
512 pharris 1.29 Double_t JetTools::dR2Mean(const PFJet *iJet,int iPFType) {
513     double lDR2Mean = 0;
514     double lSumPt2 = 0;
515     for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
516     const PFCandidate *pCand = iJet->PFCand(i0);
517     if(iPFType != -1 && pCand->PFType() != iPFType) continue;
518     lSumPt2 += pCand->Pt()*pCand->Pt();
519     double pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
520     lDR2Mean += pDR*pDR*(pCand->Pt()*pCand->Pt());
521     }
522     lDR2Mean/=lSumPt2;
523     return lDR2Mean;
524     }
525 pharris 1.18 Double_t JetTools::frac(const PFJet *iJet,Double_t iDRMax,Double_t iDRMin,Int_t iPFType) {
526     double lFrac = 0;
527     for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
528     const PFCandidate *pCand = iJet->PFCand(i0);
529     if(iPFType != -1 && pCand->PFType() != iPFType) continue;
530     Double_t pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
531     if(pDR > iDRMax) continue;
532     if(pDR < iDRMax-0.1) continue;
533 ceballos 1.25 lFrac += pCand->Pt()/iJet->RawMom().Pt();
534 pharris 1.18 }
535     return lFrac;
536     }
537     Double_t JetTools::betaStar(const PFJet *iJet,const Vertex *iVertex,const VertexCol* iVertices,Double_t iDZCut) {
538     Double_t lTotal = 0;
539     Double_t lPileup = 0;
540     for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
541     const PFCandidate* pPF = iJet->PFCand(i0);
542     const Track* pTrack = pPF->TrackerTrk();
543 pharris 1.20 //if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW
544 pharris 1.18 if(pTrack == 0) continue;
545 ceballos 1.24 lTotal += pTrack->Pt();
546 pharris 1.20 double pDZPV = fabs(pTrack->DzCorrected(*iVertex));
547     double pDZMin = pDZPV;
548 pharris 1.18 for(unsigned int i1 = 0; i1 < iVertices->GetEntries(); i1++) {
549     const Vertex *pV = iVertices->At(i1);
550 pharris 1.20 if(pV->Ndof() < 4 ||
551 ceballos 1.23 (pV->Position() - iVertex->Position()).R() < 0.02 ) continue;
552 pharris 1.20 pDZMin = TMath::Min(pDZMin,fabs(pTrack->DzCorrected(*pV)));
553 pharris 1.18 }
554 pharris 1.20 if(pDZPV > 0.2 && pDZMin < 0.2) lPileup += pTrack->Pt();
555 pharris 1.18 }
556     if(lTotal == 0) lTotal = 1;
557     return lPileup/(lTotal);
558     }
559 pharris 1.16 Bool_t JetTools::passPFLooseId(const PFJet *iJet) {
560 pharris 1.26 if(iJet->RawMom().E() == 0) return false;
561     if(iJet->NeutralHadronEnergy()/iJet->RawMom().E() > 0.99) return false;
562     if(iJet->NeutralEmEnergy()/iJet->RawMom().E() > 0.99) return false;
563 ceballos 1.27 if(iJet->NConstituents() < 2) return false;
564 pharris 1.26 if(iJet->ChargedHadronEnergy()/iJet->RawMom().E() <= 0 && fabs(iJet->Eta()) < 2.4 ) return false;
565     if(iJet->ChargedEmEnergy()/iJet->RawMom().E() > 0.99 && fabs(iJet->Eta()) < 2.4 ) return false;
566 ceballos 1.27 if(iJet->ChargedMultiplicity() < 1 && fabs(iJet->Eta()) < 2.4 ) return false;
567 pharris 1.28 if(fabs(iJet->Eta()) > 4.99) return false;
568 pharris 1.16 return true;
569     }
570    
571     /*
572     double JetTools::genFrac(const PFJet *iJet) {
573     double lTrueFrac = 0;
574     for(UInt_t i0 = 0; i0 < fParticles->GetEntries(); i0++) {
575     const MCParticle *p = fParticles->At(i0);
576     if(p->Status() != 1) continue;
577     double pDEta = iJet->Eta() - p->Eta();
578     double pDPhi = fabs(iJet->Phi()-p->Phi()); if(pDPhi > 2.*TMath::Pi() - pDPhi) pDPhi = 2.*TMath::Pi() - pDPhi;
579     double pDR = sqrt(pDEta*pDEta + pDPhi*pDPhi);
580     if(pDR > 0.5) continue;
581     lTrueFrac += p->Pt();
582     }
583     lTrueFrac/=iJet->Pt();
584     return lTrueFrac;
585     }
586     */