Reference Manual

Go to the documentation of this file.
00001 #include "RecoHI/HiJetAlgos/plugins/HiL1Subtractor.h"
00002 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00003 #include "DataFormats/JetReco/interface/PFJetCollection.h"
00004 #include "DataFormats/JetReco/interface/GenJetCollection.h"
00005 
00006 #include <vector>
00007 
00008 using namespace std;
00009 
00010 //
00011 // constants, enums and typedefs
00012 //
00013 double puCent[11] = {-5,-4,-3,-2,-1,0,1,2,3,4,5};
00014 double medianPtkt[12];
00015 
00016 //
00017 // static data member definitions
00018 //
00019 
00020 //
00021 // constructors and destructor
00022 //
00023 HiL1Subtractor::HiL1Subtractor(const edm::ParameterSet& iConfig) :
00024   src_       ( iConfig.getParameter<edm::InputTag>("src") ),
00025   jetType_       (iConfig.getParameter<std::string>("jetType") ),
00026   rhoTag_       (iConfig.getParameter<std::string>("rhoTag") )
00027 
00028 {
00029 
00030   if(jetType_ == "CaloJet"){
00031     produces<reco::CaloJetCollection >();
00032   }
00033   else if(jetType_ == "PFJet"){
00034     produces<reco::PFJetCollection >();
00035   }
00036   else if(jetType_ == "GenJet"){
00037      produces<reco::GenJetCollection >();
00038   }
00039   else{
00040     throw cms::Exception("InvalidInput") << "invalid jet type in HiL1Subtractor\n";
00041   }
00042 
00043 }
00044 
00045 
00046 HiL1Subtractor::~HiL1Subtractor()
00047 {
00048 }
00049 
00050 
00051 //
00052 // member functions
00053 //
00054 
00055 // ------------ method called to produce the data  ------------
00056 void
00057 HiL1Subtractor::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00058 {
00059 
00060   // get the input jet collection and create output jet collection
00061 
00062   // right now, identical loop for calo and PF jets, should template
00063    if(jetType_ == "GenJet"){
00064       std::auto_ptr<reco::GenJetCollection> jets( new reco::GenJetCollection);
00065       edm::Handle< edm::View<reco::GenJet> > h_jets;
00066       iEvent.getByLabel( src_, h_jets );
00067 
00068       // Grab appropriate rho, hard coded for the moment                                                                                      
00069       edm::Handle<std::vector<double> > rs;
00070       iEvent.getByLabel(edm::InputTag(rhoTag_,"rhos"),rs);
00071       //iEvent.getByLabel(edm::InputTag("ic5CaloJets","rhos"),rs);                                                                            
00072       int rsize = rs->size();
00073 
00074       for(int j = 0; j < rsize; j++){
00075          double medianpt=rs->at(j);
00076          medianPtkt[j]=medianpt;
00077       }
00078 
00079       // loop over the jets                                                                                                                   
00080       int jetsize = h_jets->size();
00081       for(int ijet = 0; ijet < jetsize; ++ijet){
00082 
00083          reco::GenJet jet = ((*h_jets)[ijet]);
00084 
00085          double jet_eta = jet.eta();
00086          double jet_et = jet.et();
00087 
00088          //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;                                                                            
00089 
00090          if(fabs(jet_eta)<=3){
00091 
00092             double rho=-999;
00093 
00094             if (jet_eta<-2.5 && jet_eta>-3.5)rho=medianPtkt[2];
00095             if (jet_eta<-1.5 && jet_eta>-2.5)rho=medianPtkt[3];
00096             if (jet_eta<-0.5 && jet_eta>-1.5)rho=medianPtkt[4];
00097             if (jet_eta<0.5 && jet_eta>-0.5)rho=medianPtkt[5];
00098             if (jet_eta<1.5 && jet_eta>0.5)rho=medianPtkt[6];
00099             if (jet_eta<2.5 && jet_eta>1.5)rho=medianPtkt[7];
00100             if (jet_eta<3.5 && jet_eta>2.5)rho=medianPtkt[8];
00101 
00102             double jet_area = jet.jetArea();
00103 
00104             double CorrFactor =0.;
00105             if(rho*jet_area<jet_et) CorrFactor = 1.0 - rho*jet_area/jet_et;
00106             jet.scaleEnergy( CorrFactor );
00107             jet.setPileup(rho*jet_area);
00108             
00109             //std::cout<<"  correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;                                                          
00110          }
00111 
00112          //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;                                                                              
00113          jets->push_back(jet);
00114 
00115 
00116       }
00117       iEvent.put(jets);
00118 
00119    }else if(jetType_ == "CaloJet"){
00120     std::auto_ptr<reco::CaloJetCollection> jets( new reco::CaloJetCollection);
00121     edm::Handle< edm::View<reco::CaloJet> > h_jets;
00122     iEvent.getByLabel( src_, h_jets );
00123 
00124     // Grab appropriate rho, hard coded for the moment
00125     edm::Handle<std::vector<double> > rs;
00126     iEvent.getByLabel(edm::InputTag(rhoTag_,"rhos"),rs);
00127     //iEvent.getByLabel(edm::InputTag("ic5CaloJets","rhos"),rs);
00128 
00129 
00130     int rsize = rs->size();
00131 
00132     for(int j = 0; j < rsize; j++){
00133       double medianpt=rs->at(j);
00134       medianPtkt[j]=medianpt;
00135     }
00136 
00137     // loop over the jets
00138 
00139     int jetsize = h_jets->size();
00140 
00141     for(int ijet = 0; ijet < jetsize; ++ijet){
00142 
00143       reco::CaloJet jet = ((*h_jets)[ijet]);
00144 
00145       double jet_eta = jet.eta();
00146       double jet_et = jet.et();
00147 
00148       //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;
00149 
00150       if(fabs(jet_eta)<=3){
00151 
00152         double rho=-999;
00153 
00154         if (jet_eta<-2.5 && jet_eta>-3.5)rho=medianPtkt[2];
00155         if (jet_eta<-1.5 && jet_eta>-2.5)rho=medianPtkt[3];
00156         if (jet_eta<-0.5 && jet_eta>-1.5)rho=medianPtkt[4];
00157         if (jet_eta<0.5 && jet_eta>-0.5)rho=medianPtkt[5];
00158         if (jet_eta<1.5 && jet_eta>0.5)rho=medianPtkt[6];
00159         if (jet_eta<2.5 && jet_eta>1.5)rho=medianPtkt[7];
00160         if (jet_eta<3.5 && jet_eta>2.5)rho=medianPtkt[8];
00161 
00162         double jet_area = jet.jetArea();
00163 
00164         double CorrFactor =0.;
00165         if(rho*jet_area<jet_et) CorrFactor = 1.0 - rho*jet_area/jet_et;
00166         jet.scaleEnergy( CorrFactor );
00167         jet.setPileup(rho*jet_area);
00168 
00169         //std::cout<<"  correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;
00170       }
00171 
00172       //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;
00173 
00174       jets->push_back(jet);
00175 
00176 
00177     }
00178     iEvent.put(jets);
00179 
00180   }
00181   else if(jetType_ == "PFJet"){
00182     std::auto_ptr<reco::PFJetCollection> jets( new reco::PFJetCollection);
00183     edm::Handle< edm::View<reco::PFJet> > h_jets;
00184     iEvent.getByLabel( src_, h_jets );
00185 
00186     // Grab appropriate rho, hard coded for the moment
00187     edm::Handle<std::vector<double> > rs;
00188     iEvent.getByLabel(edm::InputTag(rhoTag_,"rhos"),rs);
00189     //iEvent.getByLabel(edm::InputTag("ic5CaloJets","rhos"),rs);
00190 
00191 
00192     int rsize = rs->size();
00193 
00194     for(int j = 0; j < rsize; j++){
00195       double medianpt=rs->at(j);
00196       medianPtkt[j]=medianpt;
00197     }
00198 
00199     // loop over the jets
00200 
00201     int jetsize = h_jets->size();
00202 
00203     for(int ijet = 0; ijet < jetsize; ++ijet){
00204 
00205       reco::PFJet jet = ((*h_jets)[ijet]);
00206 
00207       double jet_eta = jet.eta();
00208       double jet_et = jet.et();
00209 
00210       //std::cout<<" pre-subtracted jet_et "<<jet_et<<std::endl;
00211 
00212       if(fabs(jet_eta)<=3){
00213 
00214         double rho=-999;
00215 
00216         if (jet_eta<-2.5 && jet_eta>-3.5)rho=medianPtkt[2];
00217         if (jet_eta<-1.5 && jet_eta>-2.5)rho=medianPtkt[3];
00218         if (jet_eta<-0.5 && jet_eta>-1.5)rho=medianPtkt[4];
00219         if (jet_eta<0.5 && jet_eta>-0.5)rho=medianPtkt[5];
00220         if (jet_eta<1.5 && jet_eta>0.5)rho=medianPtkt[6];
00221         if (jet_eta<2.5 && jet_eta>1.5)rho=medianPtkt[7];
00222         if (jet_eta<3.5 && jet_eta>2.5)rho=medianPtkt[8];
00223 
00224         double jet_area = jet.jetArea();
00225 
00226         double CorrFactor =0.;
00227         if(rho*jet_area<jet_et) CorrFactor = 1.0 - rho*jet_area/jet_et;
00228         jet.scaleEnergy( CorrFactor );
00229         jet.setPileup(rho*jet_area);
00230 
00231         //std::cout<<"  correction factor "<<1.0 - rho*jet_area/jet_et<<std::endl;
00232       }
00233 
00234       //std::cout<<" subtracted jet_et "<<jet.et()<<std::endl;
00235 
00236       jets->push_back(jet);
00237 
00238 
00239     }
00240     iEvent.put(jets);
00241 
00242   }
00243 
00244 
00245 
00246 
00247 
00248 
00249 
00250 }
00251 
00252 // ------------ method called once each job just before starting event loop  ------------
00253 void
00254 HiL1Subtractor::beginJob()
00255 {
00256 }
00257 
00258 // ------------ method called once each job just after ending the event loop  ------------
00259 void
00260 HiL1Subtractor::endJob() {
00261 }
00262 
00263 //define this as a plug-in
00264 DEFINE_FWK_MODULE(HiL1Subtractor);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/RecoHI/HiJetAlgos/plugins/HiL1Subtractor.cc
