00001 // -*- C++ -*- 00002 // 00003 // Package: TruthTauDecayModeProducer 00004 // Class: TruthTauDecayModeProducer 00005 // 00015 // 00016 // Original Author: Evan K. Friis, UC Davis (friis@physics.ucdavis.edu) 00017 // Created: Thu Sep 1 06:19:05 PST 2008 00018 // $Id: TruthTauDecayModeProducer.cc,v 1.2.2.1 2008/11/30 01:02:44 friis Exp $ 00019 // 00020 // 00021 00022 00023 // system include files 00024 #include <memory> 00025 00026 // user include files 00027 #include "FWCore/Framework/interface/Frameworkfwd.h" 00028 #include "FWCore/Framework/interface/EDProducer.h" 00029 00030 #include "FWCore/Framework/interface/Event.h" 00031 #include "FWCore/Framework/interface/MakerMacros.h" 00032 #include "DataFormats/Candidate/interface/Candidate.h" 00033 #include "DataFormats/TauReco/interface/PFTauDecayMode.h" 00034 #include "RecoTauTag/TauTagTools/interface/GeneratorTau.h" 00035 #include "DataFormats/Candidate/interface/Candidate.h" 00036 #include "DataFormats/JetReco/interface/GenJet.h" 00037 #include "DataFormats/HepMCCandidate/interface/GenParticle.h" 00038 #include "PhysicsTools/CandUtils/interface/AddFourMomenta.h" 00039 00040 #include <vector> 00041 00042 #include "FWCore/ParameterSet/interface/ParameterSet.h" 00043 00044 class TruthTauDecayModeProducer : public edm::EDProducer { 00045 public: 00046 struct tauObjectsHolder { 00047 std::vector<const reco::Candidate*> chargedObjects; 00048 std::vector<const reco::Candidate*> neutralObjects; 00049 }; 00050 explicit TruthTauDecayModeProducer(const edm::ParameterSet&); 00051 ~TruthTauDecayModeProducer(); 00052 00053 private: 00054 virtual void beginJob(const edm::EventSetup&) ; 00055 virtual void produce(edm::Event&, const edm::EventSetup&); 00056 virtual void endJob() ; 00057 00058 //for signal, the module takes input from a PdgIdAndStatusCandViewSelector 00059 //for background, the module takes input from a collection of GenJets 00060 bool iAmSignal_; 00061 edm::InputTag inputTag_; 00062 double leadTrackPtCut_; 00063 double leadTrackEtaCut_; 00064 double totalPtCut_; 00065 double totalEtaCut_; 00066 AddFourMomenta addP4; 00067 }; 00068 00069 TruthTauDecayModeProducer::TruthTauDecayModeProducer(const edm::ParameterSet& iConfig) 00070 { 00071 edm::LogInfo("TruthTauDecayModeProducer") << "Initializing ctor of TruthTauDecayModeProducer"; 00072 iAmSignal_ = iConfig.getParameter<bool>("iAmSignal"); 00073 inputTag_ = iConfig.getParameter<edm::InputTag>("inputTag"); 00074 leadTrackPtCut_ = iConfig.getParameter<double>("leadTrackPtCut"); 00075 leadTrackEtaCut_ = iConfig.getParameter<double>("leadTrackEtaCut"); 00076 totalPtCut_ = iConfig.getParameter<double>("totalPtCut"); 00077 totalEtaCut_ = iConfig.getParameter<double>("totalEtaCut"); 00078 //register your products 00079 edm::LogInfo("TruthTauDecayModeProducer") << "Registering products"; 00080 produces<std::vector<reco::PFTauDecayMode> >(); 00081 edm::LogInfo("TruthTauDecayModeProducer") << "TruthTauDecayModeProducer initialized"; 00082 } 00083 00084 00085 TruthTauDecayModeProducer::~TruthTauDecayModeProducer() 00086 { 00087 } 00088 00089 void 00090 TruthTauDecayModeProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) 00091 { 00092 using namespace edm; 00093 using namespace std; 00094 using namespace reco; 00095 00096 vector<tauObjectsHolder> tausToAdd_; 00097 00098 /* ********************************************** 00099 * ********** True Tau Case *********** 00100 * ********************************************** */ 00101 if (iAmSignal_) 00102 { 00103 Handle<RefToBaseVector<reco::Candidate> > decayedMCTaus; 00104 iEvent.getByLabel(inputTag_, decayedMCTaus); 00105 for(edm::RefToBaseVector<reco::Candidate>::const_iterator iterGen = decayedMCTaus->begin(); 00106 iterGen != decayedMCTaus->end(); 00107 ++iterGen) 00108 { 00109 //copy into custom format (this casting is bullshit) 00110 GeneratorTau tempTau = static_cast<const GenParticle&>(*(*iterGen)); 00111 //LogInfo("MCTauCandidateProducer") << "Generator tau produced, initializing.. "; 00112 tempTau.init(); 00113 //LogInfo("MCTauCandidateProducer") << "GenTau initialization done"; 00114 if (tempTau.isFinalStateTau()) 00115 { 00116 // Build a Tau Candidate from the information contained in the parsing class 00117 tauObjectsHolder tempTauHolder; 00118 tempTauHolder.chargedObjects = tempTau.getGenChargedPions(); 00119 tempTauHolder.neutralObjects = tempTau.getGenNeutralPions(); 00120 tausToAdd_.push_back(tempTauHolder); 00121 } 00122 } 00123 } else 00124 { 00125 /* ********************************************** 00126 * ********** QCD Case *********** 00127 * ********************************************** */ 00128 Handle<GenJetCollection> genJets; 00129 iEvent.getByLabel(inputTag_, genJets); 00130 for(GenJetCollection::const_iterator aGenJet = genJets->begin(); aGenJet != genJets->end(); ++aGenJet) 00131 { 00132 // get all constituents 00133 vector<const GenParticle*> theJetConstituents = aGenJet->getGenConstituents(); 00134 00135 tauObjectsHolder tempTauHolder; 00136 // filter the constituents 00137 for( vector<const GenParticle*>::const_iterator aCandidate = theJetConstituents.begin(); 00138 aCandidate != theJetConstituents.end(); 00139 ++aCandidate) 00140 { 00141 int pdgId = abs((*aCandidate)->pdgId()); 00142 const Candidate* theCandidate = static_cast<const Candidate*>(*aCandidate); 00143 //filter nus 00144 if (pdgId == 16 || pdgId == 12 || pdgId == 14) 00145 { 00146 //do nothing 00147 } else 00148 { 00149 // call everything charged a pion 00150 // call everything neutral a neutral pion 00151 if (theCandidate->charge() != 0) 00152 tempTauHolder.chargedObjects.push_back(theCandidate); 00153 else 00154 tempTauHolder.neutralObjects.push_back(theCandidate); 00155 } 00156 } 00157 tausToAdd_.push_back(tempTauHolder); 00158 } 00159 } 00160 00161 //output collection 00162 std::auto_ptr<vector<PFTauDecayMode> > pOut( new vector<PFTauDecayMode> ); 00163 for(vector<tauObjectsHolder>::const_iterator iTempTau = tausToAdd_.begin(); 00164 iTempTau != tausToAdd_.end(); 00165 ++iTempTau) 00166 { 00167 double leadTrackPt = 0.; 00168 double leadTrackEta = 0.; 00169 VertexCompositeCandidate chargedObjectsToAdd; 00170 const vector<const Candidate*>* chargedObjects = &(iTempTau->chargedObjects); 00171 for(vector<const Candidate*>::const_iterator iCharged = chargedObjects->begin(); 00172 iCharged != chargedObjects->end(); 00173 ++iCharged) 00174 { 00175 chargedObjectsToAdd.addDaughter(**iCharged); 00176 double trackPt = (*iCharged)->pt(); 00177 if (trackPt > leadTrackPt) 00178 { 00179 leadTrackPt = trackPt; 00180 leadTrackEta = (*iCharged)->eta(); 00181 } 00182 } 00183 //update the composite four vector 00184 addP4.set(chargedObjectsToAdd); 00185 00186 CompositeCandidate neutralPionsToAdd; 00187 const vector<const Candidate*>* neutralObjects = &(iTempTau->neutralObjects); 00188 for(vector<const Candidate*>::const_iterator iNeutral = neutralObjects->begin(); 00189 iNeutral != neutralObjects->end(); 00190 ++iNeutral) 00191 { 00192 neutralPionsToAdd.addDaughter(**iNeutral); 00193 } 00194 addP4.set(neutralPionsToAdd); 00195 00196 Particle::LorentzVector myFourVector = chargedObjectsToAdd.p4(); 00197 myFourVector += neutralPionsToAdd.p4(); 00198 00199 if(leadTrackPt > leadTrackPtCut_ && abs(leadTrackEta) < leadTrackEtaCut_ && myFourVector.pt() > totalPtCut_ && abs(myFourVector.eta()) < totalEtaCut_) 00200 { 00201 //TODO: add vertex fitting 00202 CompositeCandidate theOutliers; 00203 PFTauRef noPFTau; //empty REF to PFTau 00204 PFTauDecayMode decayModeToAdd(chargedObjectsToAdd, neutralPionsToAdd, theOutliers); 00205 decayModeToAdd.setPFTauRef(noPFTau); 00206 pOut->push_back(decayModeToAdd); 00207 } 00208 } 00209 iEvent.put(pOut); 00210 } 00211 00212 void 00213 TruthTauDecayModeProducer::beginJob(const edm::EventSetup&) 00214 { 00215 } 00216 00217 // ------------ method called once each job just after ending the event loop ------------ 00218 void 00219 TruthTauDecayModeProducer::endJob() { 00220 } 00221 00222 //define this as a plug-in 00223 DEFINE_FWK_MODULE(TruthTauDecayModeProducer);