BTauReco/src/TauMassTagInfo.cc

0001 #include "DataFormats/BTauReco/interface/TauMassTagInfo.h"
0002 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0003 #include "DataFormats/TrackReco/interface/Track.h"
0004 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0005 #include "DataFormats/Common/interface/RefVector.h"
0006
0007 using namespace edm;
0008 using namespace reco;
0009 using namespace std;
0010
0011 float reco::TauMassTagInfo::discriminator(
0012     double matching_cone, double leading_trk_pt, double signal_cone, double cluster_track_cone, double m_cut) const {
0013   float discriminator = 0.0;
0014   double invariantMass = getInvariantMass(matching_cone, leading_trk_pt, signal_cone, cluster_track_cone);
0015   if (invariantMass >= 0.0 && invariantMass < m_cut)
0016     discriminator = 1.0;
0017   return discriminator;
0018 }
0019 //
0020 // -- Set IsolatedTauTag
0021 //
0022 void reco::TauMassTagInfo::setIsolatedTauTag(const IsolatedTauTagInfoRef isolationRef) { isolatedTau = isolationRef; }
0023 //
0024 // -- Get IsolatedTauTag
0025 //
0026 const IsolatedTauTagInfoRef& reco::TauMassTagInfo::getIsolatedTauTag() const { return isolatedTau; }
0027 //
0028 // -- Set Cluster Collection
0029 //
0030 void reco::TauMassTagInfo::storeClusterTrackCollection(reco::BasicClusterRef clusterRef, float dr) {
0031   clusterMap.insert(clusterRef, dr);
0032 }
0033 //
0034 // -- Calculate 4 momentum vector from tracks only
0035 //
0036 bool reco::TauMassTagInfo::calculateTrkP4(double matching_cone,
0037                                           double leading_trk_pt,
0038                                           double signal_cone,
0039                                           math::XYZTLorentzVector& p4) const {
0040   const TrackRef leadTk = isolatedTau->leadingSignalTrack(matching_cone, leading_trk_pt);
0041   if (!leadTk) {
0042     std::cout << " TauMassTagInfo::  No Leading Track !!  " << std::endl;
0043     return false;
0044   }
0045   math::XYZVector momentum = (*leadTk).momentum();
0046   const RefVector<TrackCollection> signalTracks = isolatedTau->tracksInCone(momentum, signal_cone, 1.0);
0047   //  if (signalTracks.size() == 0 || signalTracks.size()%2 == 0) return false;
0048   if (signalTracks.empty())
0049     return false;
0050
0051   double px_inv = 0.0;
0052   double py_inv = 0.0;
0053   double pz_inv = 0.0;
0054   double e_inv = 0.0;
0055   for (RefVector<TrackCollection>::const_iterator itrack = signalTracks.begin(); itrack != signalTracks.end();
0056        itrack++) {
0057     double p = (*itrack)->p();
0058     double energy = sqrt(p * p + 0.139 * 0.139);  // use correct value!
0059     px_inv += (*itrack)->px();
0060     py_inv += (*itrack)->py();
0061     pz_inv += (*itrack)->pz();
0062     e_inv += energy;
0063   }
0064
0065   p4.SetPx(px_inv);
0066   p4.SetPy(py_inv);
0067   p4.SetPz(pz_inv);
0068   p4.SetE(e_inv);
0069
0070   return true;
0071 }
0072 //
0073 // -- Get Invariant Mass
0074 //
0075 double reco::TauMassTagInfo::getInvariantMassTrk(double matching_cone,
0076                                                  double leading_trk_pt,
0077                                                  double signal_cone) const {
0078   math::XYZTLorentzVector totalP4;
0079   if (!calculateTrkP4(matching_cone, leading_trk_pt, signal_cone, totalP4))
0080     return -1.0;
0081   return totalP4.M();
0082 }
0083 //
0084 // -- Get Invariant Mass
0085 //
0086 double reco::TauMassTagInfo::getInvariantMass(double matching_cone,
0087                                               double leading_trk_pt,
0088                                               double signal_cone,
0089                                               double track_cone) const {
0090   math::XYZTLorentzVector totalP4;
0091   if (!calculateTrkP4(matching_cone, leading_trk_pt, signal_cone, totalP4))
0092     return -1.0;
0093
0094   // Add Clusters away from tracks
0095   for (ClusterTrackAssociationCollection::const_iterator mapIter = clusterMap.begin(); mapIter != clusterMap.end();
0096        mapIter++) {
0097     const reco::BasicClusterRef& iclus = mapIter->key;
0098     float dr = mapIter->val;
0099     if (dr > track_cone) {
0100       math::XYZVector clus3Vec(iclus->x(), iclus->y(), iclus->z());
0101       double e = iclus->energy();
0102       double theta = clus3Vec.theta();
0103       double phi = clus3Vec.phi();
0104       double px = e * sin(theta) * cos(phi);
0105       double py = e * sin(theta) * sin(phi);
0106       double pz = e * cos(theta);
0107       math::XYZTLorentzVector p4(px, py, pz, e);
0108       totalP4 += p4;
0109     }
0110   }
0111   return totalP4.M();
0112 }