CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByMVAIsolation.cc

Go to the documentation of this file.
00001 #include <vector>
00002 #include <TFile.h>
00003 #include "DataFormats/Math/interface/deltaR.h"
00004 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
00005 #include "CondFormats/EgammaObjects/interface/GBRForest.h"
00006 
00007 /* class PFRecoTauDiscriminationByMVAIsolation
00008  *
00009  * Discriminates taus based on isolation deposit rings around tau
00010  *  
00011  *
00012  * 
00013  * Authors : Matthew Chan
00014  */
00015 
00016 using namespace std;
00017 using namespace reco;
00018 
00019 namespace reco {
00020   namespace tau {
00021     namespace cone {
00022       struct IsoRings
00023       {
00024         vector<int> niso;
00025         vector<vector<float> > rings;
00026         vector<vector<float> > shapes;
00027 
00028         vector<float> getVector()
00029         {
00030           vector<float> all;
00031           all.reserve(33);
00032     
00033           for(unsigned int i = 0; i < niso.size(); i++)
00034             all.push_back(niso[i]);
00035 
00036           for(unsigned int i = 0; i < rings.size(); i++)
00037             all.insert(all.end(), rings[i].begin(), rings[i].end());
00038 
00039           for(unsigned int i = 0; i < shapes.size(); i++)
00040             all.insert(all.end(), shapes[i].begin(), shapes[i].end());
00041 
00042           return all;
00043         }
00044       };
00045     }
00046   }
00047 }
00048 
00049 class PFRecoTauDiscriminationByMVAIsolation : public PFTauDiscriminationProducerBase
00050 {
00051 public:
00052   explicit PFRecoTauDiscriminationByMVAIsolation(const edm::ParameterSet& iConfig) : 
00053     PFTauDiscriminationProducerBase(iConfig),
00054     rhoProducer_(iConfig.getParameter<edm::InputTag>("rhoProducer")),
00055     gbrfFilePath_(iConfig.getParameter<edm::FileInPath>("gbrfFilePath")),
00056     returnMVA_(iConfig.getParameter<bool>("returnMVA")),
00057     mvaMin_(iConfig.getParameter<double>("mvaMin")),
00058     rho_(0)
00059   {
00060     // Prediscriminant fail value
00061     if(returnMVA_)
00062       prediscriminantFailValue_ = -1;
00063     else
00064       prediscriminantFailValue_ = 0;
00065 
00066     // Read GBRForest
00067     TFile *gbrfFile = new TFile(gbrfFilePath_.fullPath().data());
00068     gbrfTauIso_ = (GBRForest *)(gbrfFile->Get("gbrfTauIso"));
00069   }
00070 
00071   ~PFRecoTauDiscriminationByMVAIsolation(){} 
00072 
00073   void beginEvent(const edm::Event& evt, const edm::EventSetup& evtSetup);
00074   double discriminate(const PFTauRef& pfTau);
00075   reco::tau::cone::IsoRings computeIsoRings(const PFTauRef& pfTau);
00076 
00077 private:
00078   edm::InputTag rhoProducer_;
00079   edm::FileInPath gbrfFilePath_;
00080   GBRForest *gbrfTauIso_;
00081   bool returnMVA_;
00082   double mvaMin_;
00083   double rho_;
00084 };
00085 
00086 void PFRecoTauDiscriminationByMVAIsolation::beginEvent(const edm::Event& event,
00087     const edm::EventSetup& eventSetup)
00088 {
00089   // Get rho of event
00090   edm::Handle<double> hRho;
00091   event.getByLabel(rhoProducer_, hRho);
00092   rho_ = *hRho;
00093 }
00094 
00095 double PFRecoTauDiscriminationByMVAIsolation::discriminate(const PFTauRef& thePFTauRef)
00096 {
00097   reco::tau::cone::IsoRings isoRings = computeIsoRings(thePFTauRef);
00098   vector<float> mvainput = isoRings.getVector();
00099   mvainput.push_back(rho_);
00100   double mvaValue = gbrfTauIso_->GetClassifier(&mvainput[0]);
00101 
00102   return returnMVA_ ? mvaValue : mvaValue > mvaMin_;
00103 }
00104 
00105 reco::tau::cone::IsoRings PFRecoTauDiscriminationByMVAIsolation::computeIsoRings(const PFTauRef& pfTau)
00106 {
00107   vector<int>            niso(3);
00108   vector<vector<float> > rings(3, vector<float>(5));
00109   vector<vector<float> > shapes(3, vector<float>(5));
00110   vector<float>          isoptsum(3);
00111 
00112   for(unsigned int i = 0; i < pfTau->isolationPFCands().size(); i++)
00113   {
00114     const PFCandidateRef pf = pfTau->isolationPFCands().at(i);
00115 
00116     // Angular distance between PF candidate and tau
00117     float deta = pfTau->eta() - pf->eta();
00118     float dphi = reco::deltaPhi(pfTau->phi(), pf->phi());
00119     float dr = reco::deltaR(pfTau->eta(), pfTau->phi(), pf->eta(), pf->phi());
00120     int pftype = 0;
00121 
00122     // Determine PF candidate type
00123     if(pf->charge() != 0)                           pftype = 0;
00124     else if(pf->particleId() == PFCandidate::gamma) pftype = 1;
00125     else                                            pftype = 2;
00126 
00127     // Number of isolation candidates by type
00128     niso[pftype]++;
00129 
00130     // Isolation Rings
00131     if(dr < 0.1)      rings[pftype][0] += pf->pt();
00132     else if(dr < 0.2) rings[pftype][1] += pf->pt();
00133     else if(dr < 0.3) rings[pftype][2] += pf->pt();
00134     else if(dr < 0.4) rings[pftype][3] += pf->pt();
00135     else if(dr < 0.5) rings[pftype][4] += pf->pt();
00136 
00137     // Angle Shape Variables
00138     shapes[pftype][0] += pf->pt() * deta;
00139     shapes[pftype][1] += pf->pt() * dphi;
00140     shapes[pftype][2] += pf->pt() * deta*deta;
00141     shapes[pftype][3] += pf->pt() * dphi*dphi;
00142     shapes[pftype][4] += pf->pt() * deta*dphi;
00143     isoptsum[pftype]  += pf->pt();
00144   }
00145 
00146   // Mean and variance of angle variables are weighted by pT
00147   for(unsigned int i = 0; i < shapes.size(); i++)
00148   {
00149     for(unsigned int j = 0; j < shapes[i].size(); j++)
00150     {
00151       shapes[i][j] = isoptsum[i] > 0 ? fabs(shapes[i][j]/isoptsum[i]) : 0;
00152     }
00153   }
00154 
00155   // Fill IsoRing object
00156   reco::tau::cone::IsoRings isoRings;
00157   isoRings.niso = niso;
00158   isoRings.rings = rings;
00159   isoRings.shapes = shapes;
00160 
00161   return isoRings;
00162 }
00163 DEFINE_FWK_MODULE(PFRecoTauDiscriminationByMVAIsolation);