Reference Manual

Go to the documentation of this file.
00001 /*
00002  * RecoTauPiZeroStripPlugin2
00003  *
00004  * Merges PFGammas in a PFJet into Candidate piZeros defined as
00005  * strips in eta-phi.
00006  *
00007  * Author: Michail Bachtis (University of Wisconsin)
00008  *
00009  * Code modifications: Evan Friis (UC Davis),
00010  *                     Christian Veelken (LLR)
00011  *
00012  * $Id $
00013  */
00014 #include <algorithm>
00015 #include <memory>
00016 
00017 #include "RecoTauTag/RecoTau/interface/RecoTauPiZeroPlugins.h"
00018 
00019 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
00020 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
00021 #include "DataFormats/VertexReco/interface/Vertex.h"
00022 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00023 #include "DataFormats/TauReco/interface/RecoTauPiZero.h"
00024 #include "DataFormats/JetReco/interface/PFJet.h"
00025 #include "CommonTools/CandUtils/interface/AddFourMomenta.h"
00026 #include "DataFormats/Math/interface/deltaPhi.h"
00027 
00028 #include "RecoTauTag/RecoTau/interface/RecoTauCommonUtilities.h"
00029 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
00030 #include "RecoTauTag/RecoTau/interface/RecoTauVertexAssociator.h"
00031 #include "RecoTauTag/RecoTau/interface/CombinatoricGenerator.h"
00032 
00033 //-------------------------------------------------------------------------------
00034 // CV: the following headers are needed only for debug print-out
00035 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00036 #include "DataFormats/TrackReco/interface/Track.h"
00037 //-------------------------------------------------------------------------------
00038 
00039 namespace reco { namespace tau {
00040 
00041 namespace {
00042   // Apply a hypothesis on the mass of the strips.
00043   math::XYZTLorentzVector applyMassConstraint(
00044       const math::XYZTLorentzVector& vec,double mass) {
00045     double factor = sqrt(vec.energy()*vec.energy()-mass*mass)/vec.P();
00046     return math::XYZTLorentzVector(
00047         vec.px()*factor,vec.py()*factor,vec.pz()*factor,vec.energy());
00048   }
00049 }
00050 
00051 class RecoTauPiZeroStripPlugin2 : public RecoTauPiZeroBuilderPlugin 
00052 {
00053  public:
00054   explicit RecoTauPiZeroStripPlugin2(const edm::ParameterSet&);
00055   virtual ~RecoTauPiZeroStripPlugin2();
00056   // Return type is auto_ptr<PiZeroVector>
00057   return_type operator()(const reco::PFJet&) const;
00058   // Hook to update PV information
00059   virtual void beginEvent();
00060   
00061  private:
00062   typedef std::vector<reco::PFCandidatePtr> PFCandPtrs;
00063   void addCandsToStrip(RecoTauPiZero&, PFCandPtrs&, const std::vector<bool>&, std::set<size_t>&, bool&) const;
00064 
00065   RecoTauVertexAssociator vertexAssociator_;
00066 
00067   RecoTauQualityCuts* qcuts_;
00068   bool applyElecTrackQcuts_;
00069   double minGammaEtStripSeed_;
00070   double minGammaEtStripAdd_;
00071 
00072   double minStripEt_;
00073 
00074   std::vector<int> inputPdgIds_;  // type of candidates to clusterize
00075   double etaAssociationDistance_; // size of strip clustering window in eta direction
00076   double phiAssociationDistance_; // size of strip clustering window in phi direction
00077 
00078   bool updateStripAfterEachDaughter_;
00079   int maxStripBuildIterations_;
00080 
00081   // Parameters for build strip combinations
00082   bool combineStrips_;
00083   int maxStrips_;
00084   double combinatoricStripMassHypo_;
00085 
00086   AddFourMomenta p4Builder_;
00087 };
00088 
00089 RecoTauPiZeroStripPlugin2::RecoTauPiZeroStripPlugin2(const edm::ParameterSet& pset)
00090   : RecoTauPiZeroBuilderPlugin(pset),
00091     vertexAssociator_(pset.getParameter<edm::ParameterSet>("qualityCuts")),
00092     qcuts_(0)
00093 {
00094   //std::cout << "<RecoTauPiZeroStripPlugin2::RecoTauPiZeroStripPlugin2>:" << std::endl;
00095 
00096   minGammaEtStripSeed_ = pset.getParameter<double>("minGammaEtStripSeed");
00097   //std::cout << " minGammaEtStripSeed = " << minGammaEtStripSeed_ << std::endl;
00098   minGammaEtStripAdd_ = pset.getParameter<double>("minGammaEtStripAdd");
00099   //std::cout << " minGammaEtStripAdd = " << minGammaEtStripAdd_ << std::endl;
00100 
00101   minStripEt_ = pset.getParameter<double>("minStripEt");
00102   //std::cout << " minStripEt = " << minStripEt_ << std::endl;
00103   
00104   edm::ParameterSet qcuts_pset = pset.getParameterSet("qualityCuts").getParameterSet("signalQualityCuts");
00105 //-------------------------------------------------------------------------------
00106 // CV: disable track quality cuts for PFElectronsPFElectron
00107 //       (treat PFElectrons like PFGammas for the purpose of building eta-phi strips)
00108   applyElecTrackQcuts_ = pset.getParameter<bool>("applyElecTrackQcuts");
00109   if ( !applyElecTrackQcuts_ ) {
00110     qcuts_pset.addParameter<double>("minTrackPt", std::min(minGammaEtStripSeed_, minGammaEtStripAdd_));
00111     qcuts_pset.addParameter<double>("maxTrackChi2", 1.e+9);
00112     qcuts_pset.addParameter<double>("maxTransverseImpactParameter", 1.e+9);
00113     qcuts_pset.addParameter<double>("maxDeltaZ", 1.);
00114     qcuts_pset.addParameter<double>("minTrackVertexWeight", -1.);
00115     qcuts_pset.addParameter<unsigned>("minTrackPixelHits", 0);
00116     qcuts_pset.addParameter<unsigned>("minTrackHits", 0);
00117   }
00118 //-------------------------------------------------------------------------------
00119   qcuts_pset.addParameter<double>("minGammaEt", std::min(minGammaEtStripSeed_, minGammaEtStripAdd_));
00120   qcuts_ = new RecoTauQualityCuts(qcuts_pset);
00121 
00122   inputPdgIds_ = pset.getParameter<std::vector<int> >("stripCandidatesParticleIds");
00123   etaAssociationDistance_ = pset.getParameter<double>("stripEtaAssociationDistance");
00124   phiAssociationDistance_ = pset.getParameter<double>("stripPhiAssociationDistance");
00125 
00126   updateStripAfterEachDaughter_ = pset.getParameter<bool>("updateStripAfterEachDaughter");
00127   maxStripBuildIterations_ = pset.getParameter<int>("maxStripBuildIterations");
00128 
00129   combineStrips_ = pset.getParameter<bool>("makeCombinatoricStrips");
00130   if ( combineStrips_ ) {
00131     maxStrips_ = pset.getParameter<int>("maxInputStrips");
00132     combinatoricStripMassHypo_ = pset.getParameter<double>("stripMassWhenCombining");
00133   }
00134 }
00135   
00136 RecoTauPiZeroStripPlugin2::~RecoTauPiZeroStripPlugin2()
00137 {
00138   delete qcuts_;
00139 }
00140 
00141 // Update the primary vertex
00142 void RecoTauPiZeroStripPlugin2::beginEvent() 
00143 {
00144   vertexAssociator_.setEvent(*evt());
00145 }
00146 
00147 void RecoTauPiZeroStripPlugin2::addCandsToStrip(RecoTauPiZero& strip, PFCandPtrs& cands, const std::vector<bool>& candFlags, 
00148                                                 std::set<size_t>& candIdsCurrentStrip, bool& isCandAdded) const
00149 {
00150   size_t numCands = cands.size();
00151   for ( size_t candId = 0; candId < numCands; ++candId ) {
00152     if ( (!candFlags[candId]) && candIdsCurrentStrip.find(candId) == candIdsCurrentStrip.end() ) { // do not include same cand twice
00153       reco::PFCandidatePtr cand = cands[candId];
00154       if ( fabs(strip.eta() - cand->eta()) < etaAssociationDistance_ && // check if cand is within eta-phi window centered on strip 
00155            fabs(strip.phi() - cand->phi()) < phiAssociationDistance_ ) {
00156         //std::cout << "--> candId = " << candId << " has been added." << std::endl;
00157         strip.addDaughter(cand);
00158         if ( updateStripAfterEachDaughter_ ) p4Builder_.set(strip);
00159         isCandAdded = true;
00160         candIdsCurrentStrip.insert(candId);
00161       }
00162     }
00163   }
00164 }
00165 
00166 void markCandsInStrip(std::vector<bool>& candFlags, const std::set<size_t>& candIds)
00167 {
00168   for ( std::set<size_t>::const_iterator candId = candIds.begin();
00169         candId != candIds.end(); ++candId ) {
00170     candFlags[*candId] = true;
00171   }
00172 }
00173 
00174 namespace {
00175   const reco::TrackBaseRef getTrack(const PFCandidate& cand) 
00176   {
00177     if      ( cand.trackRef().isNonnull()    ) return reco::TrackBaseRef(cand.trackRef());
00178     else if ( cand.gsfTrackRef().isNonnull() ) return reco::TrackBaseRef(cand.gsfTrackRef());
00179     else return reco::TrackBaseRef();
00180   }
00181 }
00182 
00183 RecoTauPiZeroStripPlugin2::return_type RecoTauPiZeroStripPlugin2::operator()(const reco::PFJet& jet) const 
00184 {
00185   PiZeroVector output;
00186 
00187   // Get the candidates passing our quality cuts
00188   qcuts_->setPV(vertexAssociator_.associatedVertex(jet));
00189   PFCandPtrs candsVector = qcuts_->filterRefs(pfCandidates(jet, inputPdgIds_));
00190 
00191   // Convert to stl::list to allow fast deletions
00192   PFCandPtrs seedCands;
00193   PFCandPtrs addCands;
00194   int idx = 0;
00195   for ( PFCandPtrs::iterator cand = candsVector.begin();
00196         cand != candsVector.end(); ++cand ) {
00197     //std::cout << "PFGamma (" << idx << "): Et = " << (*cand)->et() << "," 
00198     //          << " eta = " << (*cand)->eta() << ", phi = " << (*cand)->phi(); 
00199     if ( (*cand)->et() > minGammaEtStripSeed_ ) {
00200       //std::cout << " --> assigning seedCandId = " << seedCands.size() << std::endl;
00201       const reco::TrackBaseRef candTrack = getTrack(*cand);
00202       if ( candTrack.isNonnull() ) {
00203         //std::cout << "has Track: pt = " << candTrack->pt() << " +/- " << candTrack->ptError() << "," 
00204         //          << " eta = " << candTrack->eta() << ", phi = " << candTrack->phi() << ","
00205         //          << " charge = " << candTrack->charge() << std::endl;
00206         //std::cout << " chi2 = " << candTrack->normalizedChi2() << std::endl;
00207         //std::cout << " dIP = " << std::abs(candTrack->dxy(vertexAssociator_.associatedVertex(jet)->position())) << std::endl;
00208         //std::cout << " dZ = " << std::abs(candTrack->dz(vertexAssociator_.associatedVertex(jet)->position())) << std::endl;
00209         //std::cout << " vtxAssocWeight = " << vertexAssociator_.associatedVertex(jet)->trackWeight(candTrack) << std::endl;
00210         //std::cout << " numPxlHits = " << candTrack->hitPattern().numberOfValidPixelHits() << std::endl;
00211         //std::cout << " numTrkHits = " << candTrack->hitPattern().numberOfValidHits() << std::endl;
00212       }
00213       //std::cout << "ECAL Et: calibrated = " << (*cand)->ecalEnergy()*sin((*cand)->theta()) << "," 
00214       //          << " raw = " << (*cand)->rawEcalEnergy()*sin((*cand)->theta()) << std::endl;
00215       //std::cout << "HCAL Et: calibrated = " << (*cand)->hcalEnergy()*sin((*cand)->theta()) << "," 
00216       //          << " raw = " << (*cand)->rawHcalEnergy()*sin((*cand)->theta()) << std::endl;
00217       seedCands.push_back(*cand);
00218     } else if ( (*cand)->et() > minGammaEtStripAdd_  ) {
00219       //std::cout << " --> assigning addCandId = " << addCands.size();
00220       addCands.push_back(*cand);
00221     }
00222     //std::cout << std::endl;
00223     ++idx;
00224   }
00225 
00226   std::vector<bool> seedCandFlags(seedCands.size()); // true/false: seedCand is already/not yet included in strip
00227   std::vector<bool> addCandFlags(addCands.size());   // true/false: addCand  is already/not yet included in strip
00228 
00229   std::set<size_t> seedCandIdsCurrentStrip;
00230   std::set<size_t> addCandIdsCurrentStrip;
00231 
00232   size_t idxSeed = 0;
00233   while ( idxSeed < seedCands.size() ) {
00234     //std::cout << "idxSeed = " << idxSeed << std::endl;
00235 
00236     seedCandIdsCurrentStrip.clear();
00237     addCandIdsCurrentStrip.clear();
00238 
00239     std::auto_ptr<RecoTauPiZero> strip(new RecoTauPiZero(*seedCands[idxSeed], RecoTauPiZero::kStrips));
00240     strip->addDaughter(seedCands[idxSeed]);
00241     seedCandIdsCurrentStrip.insert(idxSeed);
00242 
00243     bool isCandAdded;
00244     int stripBuildIteration = 0;
00245     do {
00246       isCandAdded = false;
00247 
00248       //std::cout << " adding seedCands to strip..." << std::endl;
00249       addCandsToStrip(*strip, seedCands, seedCandFlags, seedCandIdsCurrentStrip, isCandAdded);
00250       //std::cout << " adding addCands to strip..." << std::endl;
00251       addCandsToStrip(*strip, addCands,  addCandFlags,  addCandIdsCurrentStrip, isCandAdded);
00252 
00253       if ( !updateStripAfterEachDaughter_ ) p4Builder_.set(*strip);
00254 
00255       ++stripBuildIteration;
00256     } while ( isCandAdded && (stripBuildIteration < maxStripBuildIterations_ || maxStripBuildIterations_ == -1) );
00257 
00258     if ( strip->et() > minStripEt_ ) { // strip passed Et cuts, add it to the event
00259       //std::cout << "Strip: Et = " << strip->et() << "," 
00260       //          << " eta = " << strip->eta() << ", phi = " << strip->phi() 
00261       //          << " --> building it !!" << std::endl;
00262 
00263       // Update the vertex
00264       if ( strip->daughterPtr(0).isNonnull() ) strip->setVertex(strip->daughterPtr(0)->vertex());
00265       output.push_back(strip);
00266 
00267       // Mark daughters as being part of this strip
00268       markCandsInStrip(seedCandFlags, seedCandIdsCurrentStrip);
00269       markCandsInStrip(addCandFlags,  addCandIdsCurrentStrip);
00270     } else { // strip failed Et cuts, just skip it
00271       //std::cout << "Strip: Et = " << strip->et() << "," 
00272       //          << " eta = " << strip->eta() << ", phi = " << strip->phi() 
00273       //          << " --> discarding it !!" << std::endl;
00274     }
00275 
00276     ++idxSeed;
00277     while ( idxSeed < seedCands.size() && seedCandFlags[idxSeed] ) {
00278       ++idxSeed; // fast-forward to next seed cand not yet included in any strip
00279     }
00280   }
00281 
00282   // Check if we want to combine our strips
00283   if ( combineStrips_ && output.size() > 1 ) {
00284     PiZeroVector stripCombinations;
00285     // Sort the output by descending pt
00286     output.sort(output.begin(), output.end(),
00287         boost::bind(&RecoTauPiZero::pt, _1) >
00288         boost::bind(&RecoTauPiZero::pt, _2));
00289     // Get the end of interesting set of strips to try and combine
00290     PiZeroVector::const_iterator end_iter = takeNElements(
00291         output.begin(), output.end(), maxStrips_);
00292 
00293     // Look at all the combinations
00294     for ( PiZeroVector::const_iterator first = output.begin();
00295           first != end_iter-1; ++first ) {
00296       for ( PiZeroVector::const_iterator second = first+1;
00297             second != end_iter; ++second ) {
00298         Candidate::LorentzVector firstP4 = first->p4();
00299         Candidate::LorentzVector secondP4 = second->p4();
00300         // If we assume a certain mass for each strip apply it here.
00301         firstP4 = applyMassConstraint(firstP4, combinatoricStripMassHypo_);
00302         secondP4 = applyMassConstraint(secondP4, combinatoricStripMassHypo_);
00303         Candidate::LorentzVector totalP4 = firstP4 + secondP4;
00304         // Make our new combined strip
00305         std::auto_ptr<RecoTauPiZero> combinedStrips(
00306             new RecoTauPiZero(0, totalP4,
00307               Candidate::Point(0, 0, 0),
00308               //111, 10001, true, RecoTauPiZero::kCombinatoricStrips));
00309               111, 10001, true, RecoTauPiZero::kUndefined));
00310 
00311         // Now loop over the strip members
00312         BOOST_FOREACH(const RecoTauPiZero::daughters::value_type& gamma,
00313             first->daughterPtrVector()) {
00314           combinedStrips->addDaughter(gamma);
00315         }
00316         BOOST_FOREACH(const RecoTauPiZero::daughters::value_type& gamma,
00317             second->daughterPtrVector()) {
00318           combinedStrips->addDaughter(gamma);
00319         }
00320         // Update the vertex
00321         if ( combinedStrips->daughterPtr(0).isNonnull() )
00322           combinedStrips->setVertex(combinedStrips->daughterPtr(0)->vertex());
00323         // Add to our collection of combined strips
00324         stripCombinations.push_back(combinedStrips);
00325       }
00326     }
00327     // When done doing all the combinations, add the combined strips to the
00328     // output.
00329     output.transfer(output.end(), stripCombinations);
00330   }
00331 
00332   return output.release();
00333 }
00334 }} // end namespace reco::tau
00335 
00336 #include "FWCore/Framework/interface/MakerMacros.h"
00337 DEFINE_EDM_PLUGIN(RecoTauPiZeroBuilderPluginFactory,
00338     reco::tau::RecoTauPiZeroStripPlugin2, "RecoTauPiZeroStripPlugin2");
/data/refman/pasoursint/CMSSW_5_3_0/src/RecoTauTag/RecoTau/plugins/RecoTauPiZeroStripPlugin2.cc
