ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetTools.cc
(Generate patch)

Comparing UserCode/MitPhysics/Utils/src/JetTools.cc (file contents):
Revision 1.11 by ceballos, Fri Mar 11 15:13:09 2011 UTC vs.
Revision 1.28 by pharris, Tue May 1 16:53:38 2012 UTC

# Line 1 | Line 1
1   #include "MitPhysics/Utils/interface/JetTools.h"
2 + #include <algorithm>
3 + #include <vector>
4  
5   ClassImp(mithep::JetTools)
6  
7   using namespace mithep;
8 +
9  
10   JetTools::JetTools()
11   {
# Line 27 | Line 30 | Double_t JetTools::NJettiness(const Part
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 <                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi()))));
33 >                  - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi())))));
34      }
35      fval = fval + fvalpart;
36    }
# Line 49 | Line 52 | Double_t JetTools::NJettiness(const PFCa
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 <                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi()))));
55 >                  - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi())))));
56      }
57      fval = fval + fvalpart;
58    }
# Line 71 | Line 74 | Double_t JetTools::NJettiness(const Trac
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 <                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi()))));
77 >                  - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi())))));
78      }
79      fval = fval + fvalpart;
80    }
# Line 93 | Line 96 | Double_t JetTools::NJettiness(const JetO
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 <                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi()))));
99 >                  - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi())))));
100      }
101      fval = fval + fvalpart;
102    }
# Line 115 | Line 118 | Double_t JetTools::NJettiness(const Calo
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 <                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi()))));
121 >                  - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi())))));
122      }
123      fval = fval + fvalpart;
124    }
# Line 288 | Line 291 | Double_t JetTools::MtHiggs(const Particl
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    }
294 <  else if(nsel == 7){ // Use of M mass and mnu == 0, no fancy stuff
295 <    enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mass()*dilepton->Mass());
296 <    enenn = TMath::Sqrt(met->Pt() *met->Pt()  + 0.0*0.0);
297 <    enex  = dilepton->Px() + met->Px();
295 <    eney  = dilepton->Py() + met->Py();
296 <    mtHiggs = (enell+enenn)*(enell+enenn) - enex*enex - eney*eney;
294 >  else if(nsel == 7){ // Use of M mass and mnu == 0
295 >    double deltaPhiDileptonMet = fabs(MathUtils::DeltaPhi(dilepton->Phi(),
296 >                                                          met->Phi()));
297 >    mtHiggs = 2.0*dilepton->Pt()*met->Pt()*(1.0 - cos(deltaPhiDileptonMet));
298    }
299  
300    if(nsel >= 0 && nsel <= 4){
# Line 337 | Line 338 | Double_t JetTools::Beta(const TrackCol *
338  
339  
340   Double_t JetTools::Beta(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){  
341 <  double Pt_jets_X = 0. ;
342 <  double Pt_jets_Y = 0. ;
343 <  double Pt_jets_X_tot = 0. ;
344 <  double Pt_jets_Y_tot = 0. ;
341 >  double Pt_jets= 0. ;
342 >  double Pt_jetsTot = 0. ;
343 >  
344 >  for(UInt_t i=0;i<jet->NPFCands();i++){
345 >    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 >      if(pDz < delta_z){
349 >        Pt_jets += jet->PFCand(i)->TrackerTrk()->Pt();
350 >      }
351 >    }
352 >  }
353 >
354 >  Double_t beta = 0.;
355 >  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_jets_X_tot += jet->PFCand(i)->BestTrk()->Px();
348 <      Pt_jets_Y_tot += jet->PFCand(i)->BestTrk()->Py();  
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_X += jet->PFCand(i)->BestTrk()->Px();
352 <        Pt_jets_Y += jet->PFCand(i)->BestTrk()->Py();
370 >        Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
371        }
372      }
373    }
374  
375 <  if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
376 <    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);
377 <
378 <  return 1.0;
375 >  Double_t beta = 1.0;
376 >  if (Pt_jetsTot > 0)
377 >    beta = Pt_jets/Pt_jetsTot;
378 >  return beta;
379   }
380  
381  
# Line 378 | Line 396 | Bool_t  JetTools::PassBetaVertexAssociat
396    return passBetaCut;
397  
398   }
399 +
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 + 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 +    if(pCand->TrackerTrk() == 0) continue;
495 +    //if(pCand->Pt() < 1.) continue; => previous iterations
496 +    if(iDZ)  lDZCorr = fabs(pCand->TrackerTrk()->DzCorrected(*iVertex));
497 +    if(!iDZ) lDZCorr = fabs(pCand->TrackerTrk()->D0Corrected(*iVertex));
498 +    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 +    lDRMean    += pDR*(pCand->Pt())/iJet->RawMom().Pt();
509 +  }
510 +  return lDRMean;
511 + }
512 + Double_t JetTools::frac(const PFJet *iJet,Double_t iDRMax,Double_t iDRMin,Int_t iPFType) {
513 +  double lFrac = 0;
514 +  for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
515 +    const PFCandidate *pCand = iJet->PFCand(i0);
516 +    if(iPFType != -1 && pCand->PFType() != iPFType) continue;
517 +    Double_t pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
518 +    if(pDR > iDRMax) continue;
519 +    if(pDR < iDRMax-0.1) continue;
520 +    lFrac += pCand->Pt()/iJet->RawMom().Pt();
521 +  }
522 +  return lFrac;
523 + }
524 + Double_t JetTools::betaStar(const PFJet *iJet,const Vertex *iVertex,const VertexCol* iVertices,Double_t iDZCut) {
525 +  Double_t lTotal = 0;  
526 +  Double_t lPileup = 0;
527 +  for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
528 +    const PFCandidate* pPF   = iJet->PFCand(i0);
529 +    const Track* pTrack      = pPF->TrackerTrk();
530 +    //if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW
531 +    if(pTrack == 0) continue;
532 +    lTotal += pTrack->Pt();
533 +    double pDZPV  = fabs(pTrack->DzCorrected(*iVertex));
534 +    double pDZMin = pDZPV;
535 +    for(unsigned int i1 = 0; i1 < iVertices->GetEntries(); i1++) {
536 +      const Vertex *pV = iVertices->At(i1);
537 +      if(pV->Ndof() < 4 ||
538 +         (pV->Position() - iVertex->Position()).R() < 0.02 ) continue;
539 +      pDZMin = TMath::Min(pDZMin,fabs(pTrack->DzCorrected(*pV)));
540 +    }
541 +    if(pDZPV > 0.2 && pDZMin < 0.2) lPileup += pTrack->Pt();
542 +  }
543 +  if(lTotal == 0) lTotal = 1;
544 +  return lPileup/(lTotal);
545 + }
546 + Bool_t  JetTools::passPFLooseId(const PFJet *iJet) {
547 +  if(iJet->RawMom().E()                              == 0)       return false;
548 +  if(iJet->NeutralHadronEnergy()/iJet->RawMom().E()  >  0.99)    return false;
549 +  if(iJet->NeutralEmEnergy()/iJet->RawMom().E()      >  0.99)    return false;
550 +  if(iJet->NConstituents()                           <  2)       return false;
551 +  if(iJet->ChargedHadronEnergy()/iJet->RawMom().E()  <= 0     && fabs(iJet->Eta()) < 2.4 ) return false;
552 +  if(iJet->ChargedEmEnergy()/iJet->RawMom().E()      >  0.99  && fabs(iJet->Eta()) < 2.4 ) return false;
553 +  if(iJet->ChargedMultiplicity()                     < 1      && fabs(iJet->Eta()) < 2.4 ) return false;
554 +  if(fabs(iJet->Eta())                               > 4.99) return false;
555 +  return true;
556 + }
557 +
558 + /*
559 + double JetTools::genFrac(const PFJet *iJet) {
560 +  double lTrueFrac = 0;
561 +  for(UInt_t i0 = 0; i0 < fParticles->GetEntries(); i0++) {
562 +    const MCParticle *p = fParticles->At(i0);
563 +    if(p->Status() != 1) continue;
564 +    double pDEta = iJet->Eta() - p->Eta();
565 +    double pDPhi = fabs(iJet->Phi()-p->Phi()); if(pDPhi > 2.*TMath::Pi() - pDPhi) pDPhi =  2.*TMath::Pi() - pDPhi;
566 +    double pDR   = sqrt(pDEta*pDEta + pDPhi*pDPhi);
567 +    if(pDR > 0.5) continue;
568 +    lTrueFrac += p->Pt();
569 +  }
570 +  lTrueFrac/=iJet->Pt();
571 +  return lTrueFrac;
572 + }
573 + */

Diff Legend

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