Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    JetPlusTracks
00004 // Class:      JetPlusTrackProducer
00005 // 
00013 //
00014 // Original Author:  Olga Kodolova,40 R-A12,+41227671273,
00015 //         Created:  Fri Feb 19 10:14:02 CET 2010
00016 // $Id: JetPlusTrackProducer.cc,v 1.4 2010/05/05 14:00:50 kodolova Exp $
00017 //
00018 //
00019 
00020 
00021 // system include files
00022 #include <memory>
00023 
00024 // user include files
00025 #include "FWCore/Framework/interface/Frameworkfwd.h"
00026 #include "FWCore/Framework/interface/EDProducer.h"
00027 
00028 #include "FWCore/Framework/interface/Event.h"
00029 #include "FWCore/Framework/interface/MakerMacros.h"
00030 
00031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00032 
00033 #include "RecoJets/JetPlusTracks/plugins/JetPlusTrackProducer.h"
00034 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00035 #include "DataFormats/JetReco/interface/CaloJet.h"
00036 #include "DataFormats/JetReco/interface/JPTJetCollection.h"
00037 #include "DataFormats/JetReco/interface/JPTJet.h"
00038 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00039 #include "DataFormats/TrackReco/interface/Track.h"
00040 #include "DataFormats/JetReco/interface/Jet.h"
00041 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00042 #include "DataFormats/VertexReco/interface/Vertex.h"
00043 #include "DataFormats/Math/interface/deltaPhi.h"
00044 #include "DataFormats/Math/interface/deltaR.h"
00045 
00046 #include <string>
00047 
00048 using namespace std;
00049 using namespace jpt;
00050 
00051 //
00052 // constants, enums and typedefs
00053 //
00054 
00055 
00056 //
00057 // static data member definitions
00058 //
00059 
00060 //
00061 // constructors and destructor
00062 //
00063 JetPlusTrackProducer::JetPlusTrackProducer(const edm::ParameterSet& iConfig)
00064 {
00065    //register your products
00066    src = iConfig.getParameter<edm::InputTag>("src");
00067    alias = iConfig.getUntrackedParameter<string>("alias");
00068    srcPVs_ = iConfig.getParameter<edm::InputTag>("srcPVs");
00069    vectorial_ = iConfig.getParameter<bool>("VectorialCorrection");
00070    useZSP = iConfig.getParameter<bool>("UseZSP");
00071    mJPTalgo  = new JetPlusTrackCorrector(iConfig);
00072    mZSPalgo  = new ZSPJPTJetCorrector(iConfig);
00073  
00074    produces<reco::JPTJetCollection>().setBranchAlias(alias); 
00075    
00076 }
00077 
00078 
00079 JetPlusTrackProducer::~JetPlusTrackProducer()
00080 {
00081  
00082    // do anything here that needs to be done at desctruction time
00083    // (e.g. close files, deallocate resources etc.)
00084 
00085 }
00086 
00087 
00088 //
00089 // member functions
00090 //
00091 
00092 // ------------ method called to produce the data  ------------
00093 void
00094 JetPlusTrackProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00095 {
00096    using namespace edm; 
00097 
00098 
00099 //  std::cout<<" RecoJets::JetPlusTrackProducer::produce "<<std::endl;
00100 
00101 
00102 // get stuff from Event
00103   edm::Handle <edm::View <reco::CaloJet> > jets_h;
00104   iEvent.getByLabel (src, jets_h);
00105 
00106 //  std::auto_ptr<reco::CaloJetCollection> pOut(new reco::CaloJetCollection());
00107   std::auto_ptr<reco::JPTJetCollection> pOut(new reco::JPTJetCollection());
00108 
00109   for (unsigned i = 0; i < jets_h->size(); ++i) {
00110 
00111    const reco::CaloJet* oldjet = &(*(jets_h->refAt(i)));
00112    
00113    reco::CaloJet corrected = *oldjet; 
00114    
00115 // ZSP corrections    
00116 
00117    double factorZSP = 1.;
00118    if(useZSP) factorZSP = mZSPalgo->correction(corrected, iEvent, iSetup);
00119 
00120    corrected.scaleEnergy (factorZSP);
00121 
00122 // JPT corrections 
00123 
00124    double scaleJPT = 1.; 
00125 
00126    math::XYZTLorentzVector p4;
00127    
00128    if ( !vectorial_ ) {
00129             
00130    scaleJPT = mJPTalgo->correction ( corrected, *oldjet, iEvent, iSetup );
00131    p4 = math::XYZTLorentzVector( corrected.px()*scaleJPT, 
00132                                  corrected.py()*scaleJPT,
00133                                  corrected.pz()*scaleJPT, 
00134                                  corrected.energy()*scaleJPT );
00135    } else {
00136    scaleJPT = mJPTalgo->correction( corrected, *oldjet, iEvent, iSetup, p4 );
00137   }         
00138 
00139 // Construct JPTJet constituent
00140   jpt::MatchedTracks pions;
00141   jpt::MatchedTracks muons;
00142   jpt::MatchedTracks elecs;
00143   
00144   bool ok = mJPTalgo->matchTracks( *oldjet, 
00145                                     iEvent, 
00146                                     iSetup,
00147                                      pions, 
00148                                      muons, 
00149                                      elecs );
00150 
00151    
00152   reco::JPTJet::Specific specific;
00153 
00154   if(ok) {
00155     specific.pionsInVertexInCalo = pions.inVertexInCalo_;
00156     specific.pionsInVertexOutCalo = pions.inVertexOutOfCalo_;
00157     specific.pionsOutVertexInCalo = pions.outOfVertexInCalo_;
00158     specific.muonsInVertexInCalo = muons.inVertexInCalo_;
00159     specific.muonsInVertexOutCalo = muons.inVertexOutOfCalo_;
00160     specific.muonsOutVertexInCalo = muons.outOfVertexInCalo_;
00161     specific.elecsInVertexInCalo = elecs.inVertexInCalo_;
00162     specific.elecsInVertexOutCalo = elecs.inVertexOutOfCalo_;
00163     specific.elecsOutVertexInCalo = elecs.outOfVertexInCalo_;
00164   }
00165 
00166 // Fill JPT Specific
00167     edm::RefToBase<reco::Jet> myjet = (edm::RefToBase<reco::Jet>)jets_h->refAt(i);
00168     specific.theCaloJetRef = myjet;
00169     specific.mZSPCor = factorZSP;
00170     specific.mResponseOfChargedWithEff = (float)mJPTalgo->getResponseOfChargedWithEff();
00171     specific.mResponseOfChargedWithoutEff = (float)mJPTalgo->getResponseOfChargedWithoutEff();
00172     specific.mSumPtOfChargedWithEff = (float)mJPTalgo->getSumPtWithEff();
00173     specific.mSumPtOfChargedWithoutEff = (float)mJPTalgo->getSumPtWithoutEff();
00174     specific.mSumEnergyOfChargedWithEff = (float)mJPTalgo->getSumEnergyWithEff();
00175     specific.mSumEnergyOfChargedWithoutEff = (float)mJPTalgo->getSumEnergyWithoutEff();
00176     specific.mChargedHadronEnergy = (float)mJPTalgo->getSumEnergyWithoutEff();
00177 
00178 // Fill Charged Jet shape parameters
00179    double deR2Tr = 0.;
00180    double deEta2Tr = 0.;
00181    double dePhi2Tr = 0.;
00182    double Zch = 0.;
00183    double Pout2 = 0.;
00184    double Pout = 0.;
00185    double denominator_tracks = 0.;
00186    int ntracks = 0;
00187 
00188    for( reco::TrackRefVector::const_iterator it = pions.inVertexInCalo_.begin(); it != pions.inVertexInCalo_.end(); it++) { 
00189     double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
00190     double deEta = (*it)->eta() - p4.eta();
00191     double dePhi = deltaPhi((*it)->phi(), p4.phi());
00192      if((**it).ptError()/(**it).pt() < 0.1) {
00193        deR2Tr   =  deR2Tr + deR*deR*(*it)->pt();
00194        deEta2Tr = deEta2Tr + deEta*deEta*(*it)->pt();
00195        dePhi2Tr = dePhi2Tr + dePhi*dePhi*(*it)->pt();
00196        denominator_tracks = denominator_tracks + (*it)->pt();
00197        Zch    =  Zch +
00198        ((*it)->px()*p4.Px()+(*it)->py()*p4.Py()+(*it)->pz()*p4.Pz())/(p4.P()*p4.P());
00199        Pout2 = Pout2 + (**it).p()*(**it).p() - (Zch*p4.P())*(Zch*p4.P());
00200        ntracks++;
00201      }
00202    }
00203 
00204 
00205 
00206    for( reco::TrackRefVector::const_iterator it = muons.inVertexInCalo_.begin(); it != muons.inVertexInCalo_.end(); it++) {
00207     double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
00208     double deEta = (*it)->eta() - p4.eta();
00209     double dePhi = deltaPhi((*it)->phi(), p4.phi());
00210      if((**it).ptError()/(**it).pt() < 0.1) {
00211        deR2Tr   =  deR2Tr + deR*deR*(*it)->pt();
00212        deEta2Tr = deEta2Tr + deEta*deEta*(*it)->pt();
00213        dePhi2Tr = dePhi2Tr + dePhi*dePhi*(*it)->pt();
00214        denominator_tracks = denominator_tracks + (*it)->pt();
00215        Zch    = Zch +
00216        ((*it)->px()*p4.Px()+(*it)->py()*p4.Py()+(*it)->pz()*p4.Pz())/(p4.P()*p4.P());
00217        Pout2 = Pout2 + (**it).p()*(**it).p() - (Zch*p4.P())*(Zch*p4.P());
00218        ntracks++;
00219      }
00220    }
00221    for( reco::TrackRefVector::const_iterator it = elecs.inVertexInCalo_.begin(); it != elecs.inVertexInCalo_.end(); it++) {
00222     double deR = deltaR((*it)->eta(), (*it)->phi(), p4.eta(), p4.phi());
00223     double deEta = (*it)->eta() - p4.eta();
00224     double dePhi = deltaPhi((*it)->phi(), p4.phi());
00225      if((**it).ptError()/(**it).pt() < 0.1) {
00226        deR2Tr   =  deR2Tr + deR*deR*(*it)->pt();
00227        deEta2Tr = deEta2Tr + deEta*deEta*(*it)->pt();
00228        dePhi2Tr = dePhi2Tr + dePhi*dePhi*(*it)->pt();
00229        denominator_tracks = denominator_tracks + (*it)->pt();
00230        Zch    = Zch + 
00231        ((*it)->px()*p4.Px()+(*it)->py()*p4.Py()+(*it)->pz()*p4.Pz())/(p4.P()*p4.P());
00232        Pout2 = Pout2 + (**it).p()*(**it).p() - (Zch*p4.P())*(Zch*p4.P());
00233        ntracks++;
00234      }
00235    }
00236 
00237         if(ntracks > 0) {
00238           Pout   = sqrt(fabs(Pout2))/ntracks;
00239           Zch    = Zch/ntracks;
00240         }
00241           if (denominator_tracks!=0){
00242             deR2Tr  = deR2Tr/denominator_tracks;
00243             deEta2Tr= deEta2Tr/denominator_tracks;
00244             dePhi2Tr= dePhi2Tr/denominator_tracks;
00245           }
00246    
00247       specific.R2momtr = deR2Tr;
00248       specific.Eta2momtr = deEta2Tr;
00249       specific.Phi2momtr = dePhi2Tr;
00250       specific.Pout = Pout;
00251       specific.Zch = Zch;
00252 
00253 
00254 //       std::cout<<" Moments for charged component "<<deR2_Tr<<" "<<deEta2_Tr<<" "<<dePhi2_Tr<<std::endl;
00255 
00256 
00257 // Create JPT jet
00258 
00259    reco::Particle::Point vertex_=reco::Jet::Point(0,0,0);
00260    
00261 // If we add primary vertex
00262    edm::Handle<reco::VertexCollection> pvCollection;
00263    iEvent.getByLabel(srcPVs_,pvCollection);
00264    if ( pvCollection.isValid() && pvCollection->size()>0 ) vertex_=pvCollection->begin()->position();
00265 
00266    reco::JPTJet fJet(p4, vertex_, specific, corrected.getJetConstituents()); 
00267 
00268   // fJet.printJet();
00269 
00270 // Output module
00271     pOut->push_back(fJet); 
00272           
00273   }
00274   
00275   iEvent.put(pOut);
00276    
00277 }
00278 
00279 // ------------ method called once each job just before starting event loop  ------------
00280 void 
00281 JetPlusTrackProducer::beginJob()
00282 {
00283 }
00284 
00285 // ------------ method called once each job just after ending the event loop  ------------
00286 void 
00287 JetPlusTrackProducer::endJob() {
00288 }
00289 
00290 //define this as a plug-in
00291 //DEFINE_FWK_MODULE(JetPlusTrackProducer);
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/RecoJets/JetPlusTracks/plugins/JetPlusTrackProducer.cc
