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.8 by ceballos, Wed Jan 19 10:18:19 2011 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 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    }
# Line 49 | Line 54 | Double_t JetTools::NJettiness(const PFCa
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(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi()))));
57 >                  - 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi())))));
58      }
59      fval = fval + fvalpart;
60    }
# Line 71 | 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    }
# Line 93 | 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    }
# Line 115 | 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    }
# Line 288 | Line 293 | Double_t JetTools::MtHiggs(const Particl
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 <= 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 300 | Line 311 | Double_t JetTools::MtHiggs(const Particl
311    return mtHiggs;
312   }
313  
314 < void JetTools::Alpha(Double_t AlphaVar[2], const TrackCol *tracks, Jet *jet, const VertexCol *vertices, Double_t  delta_z, Double_t delta_cone){  
315 <  AlphaVar[0] = -1.0;
316 <  AlphaVar[1] = -1.0;
306 <  if(tracks->GetEntries() <= 0) return;
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. ;
# Line 314 | Line 324 | void JetTools::Alpha(Double_t AlphaVar[2
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(*vertices->At(0)));
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();  
# Line 322 | Line 332 | void JetTools::Alpha(Double_t AlphaVar[2
332      }
333    }
334  
325 –  if(jet->Pt() > 0)
326 –    AlphaVar[0] = sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / jet->Pt();
335    if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
336 <    AlphaVar[1] = 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);
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 < void JetTools::Alpha(Double_t AlphaVar[2], const PFJet *jet, const VertexCol *vertices, Double_t  delta_z){  
343 <  AlphaVar[0] = -1.0;
344 <  AlphaVar[1] = -1.0;
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 Pt_jets_X = 0. ;
364 <  double Pt_jets_Y = 0. ;
365 <  double Pt_jets_X_tot = 0. ;
339 <  double Pt_jets_Y_tot = 0. ;
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_jets_X_tot += jet->PFCand(i)->BestTrk()->Px();
370 <      Pt_jets_Y_tot += jet->PFCand(i)->BestTrk()->Py();  
345 <      double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(0)));
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_X += jet->PFCand(i)->BestTrk()->Px();
348 <        Pt_jets_Y += jet->PFCand(i)->BestTrk()->Py();  
372 >        Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
373        }
374      }
375    }
376  
377 <  if(jet->Pt() > 0)
378 <    AlphaVar[0] = sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / jet->Pt();
379 <  if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
380 <    AlphaVar[1] = 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 >  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