00001 #include "RecoTauTag/InvariantMass/interface/InvariantMassAlgorithm.h"
00002 #include <boost/regex.hpp>
00003
00004 using namespace std;
00005 using namespace edm;
00006 using namespace reco;
00007
00008
00009
00010
00011 InvariantMassAlgorithm::InvariantMassAlgorithm()
00012 {
00013 matching_cone = 0.1;
00014 leading_trk_pt = 1.0;
00015 signal_cone = 0.07;
00016 cluster_jet_matching_cone = 0.4;
00017 cluster_track_matching_cone = 0.08;
00018 inv_mass_cut = 2.5;
00019 }
00020
00021
00022
00023 InvariantMassAlgorithm::InvariantMassAlgorithm(const ParameterSet& parameters)
00024 {
00025 matching_cone = parameters.getParameter<double>("MatchingCone");
00026 leading_trk_pt = parameters.getParameter<double>("LeadingTrackPt");
00027 signal_cone = parameters.getParameter<double>("SignalCone");
00028 cluster_jet_matching_cone = parameters.getParameter<double>("ClusterSelectionCone");
00029 cluster_track_matching_cone = parameters.getParameter<double>("ClusterTrackMatchingCone");
00030 inv_mass_cut = parameters.getParameter<double>("InvariantMassCutoff");
00031
00032
00033
00034 edm::ParameterSet tk_ass_pset = parameters.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
00035 trackAssociatorParameters_.loadParameters( tk_ass_pset );
00036
00037 trackAssociator_ = new TrackDetectorAssociator();
00038 trackAssociator_->useDefaultPropagator();
00039 }
00040
00041
00042
00043 InvariantMassAlgorithm::~InvariantMassAlgorithm() {
00044 if(trackAssociator_) delete trackAssociator_;
00045 }
00046
00047
00048
00049 pair<double, reco::TauMassTagInfo> InvariantMassAlgorithm::tag(edm::Event& theEvent,
00050 const edm::EventSetup& theEventSetup,
00051 const reco::IsolatedTauTagInfoRef& tauRef,
00052 const Handle<BasicClusterCollection>& clus_handle)
00053 {
00054
00055 TauMassTagInfo resultExtended;
00056 resultExtended.setIsolatedTauTag(tauRef);
00057
00058 double discriminator = tauRef->discriminator();
00059 if (discriminator > 0.0) {
00060
00061 int nSel = 0;
00062 for (size_t ic = 0; ic < clus_handle->size(); ic++) {
00063 math::XYZVector cluster3Vec((*clus_handle)[ic].x(),(*clus_handle)[ic].y(),(*clus_handle)[ic].z());
00064 BasicClusterRef cluster(clus_handle, ic);
00065 float et = (*clus_handle)[ic].energy() * sin(cluster3Vec.theta());
00066 if ( (et < 0.04) || (et > 300.0) ) continue;
00067 float delR = getMinimumClusterDR(theEvent, theEventSetup, tauRef, cluster3Vec);
00068 if (delR != -1.0) {
00069 nSel++;
00070 resultExtended.storeClusterTrackCollection(cluster, delR);
00071 }
00072 }
00073 if (0) cout << " Total Number of Clusters " << clus_handle->size() << " " << nSel << endl;
00074
00075 discriminator = resultExtended.discriminator(matching_cone,
00076 leading_trk_pt,
00077 signal_cone,
00078 cluster_track_matching_cone,
00079 inv_mass_cut);
00080 }
00081 return make_pair(discriminator, resultExtended);
00082 }
00083
00084
00085
00086 float InvariantMassAlgorithm::getMinimumClusterDR(edm::Event& theEvent,
00087 const edm::EventSetup& theEventSetup,
00088 const reco::IsolatedTauTagInfoRef& tauRef,
00089 const math::XYZVector& cluster_3vec)
00090 {
00091
00092 trackAssociatorParameters_.useEcal = true;
00093 trackAssociatorParameters_.useHcal = false;
00094 trackAssociatorParameters_.useHO = false;
00095 trackAssociatorParameters_.useMuon = false;
00096 trackAssociatorParameters_.dREcal = 0.03;
00097
00098 const TrackRefVector tracks = tauRef->allTracks();
00099 const Jet & jet = *(tauRef->jet());
00100 math::XYZVector jet3Vec(jet.px(),jet.py(),jet.pz());
00101 float min_dR = 999.9;
00102 float deltaR1 = ROOT::Math::VectorUtil::DeltaR(cluster_3vec, jet3Vec);
00103 if (deltaR1 > cluster_jet_matching_cone) return -1.0;
00104
00105 for (edm::RefVector<reco::TrackCollection>::const_iterator it = tracks.begin();
00106 it != tracks.end(); it++)
00107 {
00108 TrackDetMatchInfo info;
00109 info = trackAssociator_->associate(theEvent, theEventSetup,
00110 trackAssociator_->getFreeTrajectoryState(theEventSetup, (*(*it))),
00111 trackAssociatorParameters_);
00112 math::XYZVector track3Vec(info.trkGlobPosAtEcal.x(),
00113 info.trkGlobPosAtEcal.y(),
00114 info.trkGlobPosAtEcal.z());
00115
00116 float deltaR2 = ROOT::Math::VectorUtil::DeltaR(cluster_3vec, track3Vec);
00117 if (deltaR2 < min_dR) min_dR = deltaR2;
00118 }
00119 return min_dR;
00120 }