CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/PhysicsTools/JetMCAlgos/plugins/TauGenJetProducer.cc

Go to the documentation of this file.
00001 #include "PhysicsTools/JetMCAlgos/plugins/TauGenJetProducer.h"
00002 
00003 #include "FWCore/Framework/interface/ESHandle.h"
00004 
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 #include "FWCore/Utilities/interface/Exception.h"
00007 #include "FWCore/Framework/interface/EventSetup.h"
00008 
00009 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00010 #include "DataFormats/JetReco/interface/GenJetCollection.h"
00011 #include "DataFormats/JetReco/interface/GenJet.h"
00012 #include "DataFormats/Common/interface/RefToPtr.h"
00013 
00014 #include "PhysicsTools/HepMCCandAlgos/interface/GenParticlesHelper.h"
00015 
00016 using namespace std;
00017 using namespace edm;
00018 using namespace reco;
00019 
00020 TauGenJetProducer::TauGenJetProducer(const edm::ParameterSet& iConfig) 
00021 {
00022   inputTagGenParticles_ 
00023     = iConfig.getParameter<InputTag>("GenParticles");
00024 
00025   includeNeutrinos_ 
00026     = iConfig.getParameter<bool>("includeNeutrinos");
00027 
00028   verbose_ = 
00029     iConfig.getUntrackedParameter<bool>("verbose",false);
00030 
00031   produces<GenJetCollection>();
00032 }
00033 
00034 TauGenJetProducer::~TauGenJetProducer() { }
00035 
00036 void TauGenJetProducer::beginJob() { }
00037 
00038 void TauGenJetProducer::produce(Event& iEvent, 
00039                                 const EventSetup& iSetup) {
00040   
00041   Handle<GenParticleCollection> genParticles;
00042 
00043   bool found = iEvent.getByLabel( inputTagGenParticles_, genParticles);
00044   
00045   if ( !found ) {
00046     std::ostringstream  err;
00047     err<<" cannot get collection: "
00048        <<inputTagGenParticles_<<std::endl;
00049     edm::LogError("TauGenJetProducer")<<err.str();
00050     throw cms::Exception( "MissingProduct", err.str());
00051   }
00052 
00053   std::auto_ptr<GenJetCollection> 
00054     pOutVisTaus(new GenJetCollection());
00055 
00056   using namespace GenParticlesHelper;
00057 
00058   GenParticleRefVector allStatus2Taus;  
00059   findParticles( *genParticles,
00060                  allStatus2Taus, 15, 2);
00061 
00062   for ( IGR iTau=allStatus2Taus.begin(); iTau!=allStatus2Taus.end(); ++iTau ) {
00063 
00064     // look for all status 1 (stable) descendents 
00065     GenParticleRefVector descendents;
00066     findDescendents( *iTau, descendents, 1);
00067 
00068     // CV: skip status 2 taus that radiate-off a photon
00069     //    --> have a status 2 tau lepton in the list of descendents
00070     GenParticleRefVector status2TauDaughters;
00071     findDescendents( *iTau, status2TauDaughters, 2, 15 );
00072     if ( status2TauDaughters.size() > 0 ) continue;
00073     
00074     // loop on descendents, and take all except neutrinos
00075     math::XYZTLorentzVector sumVisMom;
00076     Particle::Charge charge = 0;
00077     Jet::Constituents constituents;
00078 
00079     if(verbose_)
00080       cout<<"tau "<<(*iTau)<<endl;
00081 
00082     for(IGR igr = descendents.begin(); 
00083         igr!= descendents.end(); ++igr ) {
00084       
00085       int absPdgId = abs((*igr)->pdgId());
00086       
00087       // neutrinos
00088       if(!includeNeutrinos_ ) {
00089         if( absPdgId == 12 || 
00090             absPdgId == 14 || 
00091             absPdgId == 16  ) 
00092           continue;
00093       }      
00094 
00095       if(verbose_) 
00096         cout<<"\t"<<(*igr)<<endl;
00097       
00098       charge += (*igr)->charge();
00099       sumVisMom += (*igr)->p4();
00100 
00101       // need to convert the vector of reference to the constituents
00102       // to a vector of pointers to build the genjet
00103       constituents.push_back( refToPtr( *igr) );
00104     }
00105 
00106     math::XYZPoint vertex;
00107     GenJet::Specific specific;
00108     
00109     GenJet jet( sumVisMom, vertex, specific, constituents); 
00110     
00111     if (charge != (*iTau)->charge() )
00112       std::cout<<" charge of Tau: " << (*iTau) << " not equal to charge of sum of charge of all descendents. " << std::cout;
00113 
00114     jet.setCharge(charge);
00115     pOutVisTaus->push_back( jet );
00116     
00117   }
00118   iEvent.put( pOutVisTaus );
00119 }
00120 
00121 
00122 DEFINE_FWK_MODULE( TauGenJetProducer );