Reference Manual

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2012/10/10 04:00:00 $
00005  *  $Revision: 1.32 $
00006  *  \author Michael B. Anderson, University of Wisconsin Madison
00007  */
00008 
00009 #include "DQM/Physics/src/QcdPhotonsDQM.h"
00010 
00011 #include "FWCore/Framework/interface/Event.h"
00012 #include "FWCore/Framework/interface/EventSetup.h"
00013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00014 
00015 #include "DQMServices/Core/interface/DQMStore.h"
00016 #include "DQMServices/Core/interface/MonitorElement.h"
00017 #include "FWCore/ServiceRegistry/interface/Service.h"
00018 #include "FWCore/Common/interface/TriggerNames.h"
00019 
00020 #include "DataFormats/Common/interface/Handle.h"
00021 
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023 
00024 // Physics Objects
00025 #include "DataFormats/EgammaCandidates/interface/Photon.h"
00026 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
00027 #include "DataFormats/JetReco/interface/Jet.h"
00028 
00029 // Vertex
00030 #include "DataFormats/VertexReco/interface/Vertex.h"
00031 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00032 
00033 // For removing ECAL Spikes
00034 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
00035 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00036 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00037 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00038 #include "FWCore/ServiceRegistry/interface/Service.h"
00039 
00040 #include "FWCore/Framework/interface/ESHandle.h"
00041 #include "CondFormats/EcalObjects/interface/EcalCondObjectContainer.h"
00042 
00043 //geometry
00044 //#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00045 //#include "Geometry/Records/interface/IdealGeometryRecord.h"
00046 //#include "Geometry/Records/interface/CaloGeometryRecord.h"
00047 //#include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00048 //#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00049 //#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00050 //#include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
00051 
00052 // Math stuff
00053 #include "DataFormats/Math/interface/deltaR.h"
00054 #include "DataFormats/Math/interface/deltaPhi.h"
00055 
00056 #include <vector>
00057 
00058 #include <string>
00059 #include <cmath>
00060 using namespace std;
00061 using namespace edm;
00062 using namespace reco;
00063 
00064 
00065 QcdPhotonsDQM::QcdPhotonsDQM(const ParameterSet& parameters) {
00066   // Get parameters from configuration file
00067   theTriggerPathToPass_        = parameters.getParameter<string>("triggerPathToPass");
00068   thePlotTheseTriggersToo_     = parameters.getParameter<vector<string> >("plotTheseTriggersToo");
00069   trigTag_                    = parameters.getUntrackedParameter<edm::InputTag>("trigTag");
00070   thePhotonCollectionLabel_    = parameters.getParameter<InputTag>("photonCollection");
00071   theJetCollectionLabel_       = parameters.getParameter<InputTag>("jetCollection");
00072   theVertexCollectionLabel_    = parameters.getParameter<InputTag>("vertexCollection");
00073   theMinJetPt_             = parameters.getParameter<double>("minJetPt");
00074   theMinPhotonEt_              = parameters.getParameter<double>("minPhotonEt");
00075   theRequirePhotonFound_       = parameters.getParameter<bool>("requirePhotonFound");
00076   thePlotPhotonMaxEt_          = parameters.getParameter<double>("plotPhotonMaxEt");
00077   thePlotPhotonMaxEta_         = parameters.getParameter<double>("plotPhotonMaxEta");
00078   thePlotJetMaxEta_            = parameters.getParameter<double>("plotJetMaxEta");
00079   theBarrelRecHitTag           = parameters.getParameter<InputTag>("barrelRecHitTag");
00080   theEndcapRecHitTag           = parameters.getParameter<InputTag>("endcapRecHitTag");
00081   // just to initialize
00082   isValidHltConfig_ = false;
00083 
00084   // coverity says...
00085   h_deltaEt_photon_jet = 0;
00086   h_deltaPhi_jet_jet2 = 0;
00087   h_deltaPhi_photon_jet = 0;
00088   h_deltaPhi_photon_jet2 = 0;
00089   h_deltaR_jet_jet2 = 0;
00090   h_deltaR_photon_jet2 = 0;
00091   h_jet2_eta = 0;
00092   h_jet2_pt = 0;
00093   h_jet2_ptOverPhotonEt = 0;
00094   h_jet_count = 0;
00095   h_jet_eta = 0;
00096   h_jet_pt = 0;
00097   h_photon_count_bar = 0;
00098   h_photon_count_end = 0;
00099   h_photon_et = 0;
00100   h_photon_et_beforeCuts = 0;
00101   h_photon_et_jetco = 0;
00102   h_photon_et_jetcs = 0;
00103   h_photon_et_jetfo = 0;
00104   h_photon_et_jetfs = 0;
00105   h_photon_et_ratio_co_cs = 0;
00106   h_photon_et_ratio_co_fo = 0;
00107   h_photon_et_ratio_co_fs = 0;
00108   h_photon_et_ratio_cs_fo = 0;
00109   h_photon_et_ratio_cs_fs = 0;
00110   h_photon_et_ratio_fo_fs = 0;
00111   h_photon_eta = 0;
00112   h_triggers_passed = 0;
00113 
00114   theDbe = Service<DQMStore>().operator->();
00115   
00116 }
00117 
00118 QcdPhotonsDQM::~QcdPhotonsDQM() { 
00119 }
00120 
00121 
00122 void QcdPhotonsDQM::beginJob() {
00123  
00124   logTraceName = "QcdPhotonAnalyzer";
00125 
00126   LogTrace(logTraceName)<<"Parameters initialization";
00127  
00128   theDbe->setCurrentFolder("Physics/QcdPhotons");  // Use folder with name of PAG
00129 
00130   std::stringstream aStringStream;
00131   std::string aString;
00132   aStringStream << theMinJetPt_;
00133   aString = aStringStream.str();
00134 
00135   // Monitor of triggers passed
00136   int numOfTriggersToMonitor = thePlotTheseTriggersToo_.size();
00137   h_triggers_passed = theDbe->book1D("triggers_passed", "Events passing these trigger paths", numOfTriggersToMonitor, 0, numOfTriggersToMonitor);
00138   for (int i=0; i<numOfTriggersToMonitor; i++) {
00139     h_triggers_passed->setBinLabel(i+1,thePlotTheseTriggersToo_[i]);
00140   }
00141 
00142   // Keep the number of plots and number of bins to a minimum!
00143   h_photon_et_beforeCuts           = theDbe->book1D("photon_et_beforeCuts",       "#gamma with highest E_{T};E_{T}(#gamma) (GeV)", 20, 0., thePlotPhotonMaxEt_);
00144   h_photon_et           = theDbe->book1D("photon_et",       "#gamma with highest E_{T};E_{T}(#gamma) (GeV)", 20, 0., thePlotPhotonMaxEt_);
00145   h_photon_eta          = theDbe->book1D("photon_eta",      "#gamma with highest E_{T};#eta(#gamma)", 40, -thePlotPhotonMaxEta_, thePlotPhotonMaxEta_);
00146   h_photon_count_bar    = theDbe->book1D("photon_count_bar","Number of #gamma's passing selection (Barrel);Number of #gamma's", 8, -0.5, 7.5);
00147   h_photon_count_end    = theDbe->book1D("photon_count_end","Number of #gamma's passing selection (Endcap);Number of #gamma's", 8, -0.5, 7.5);
00148 
00149   h_jet_pt             = theDbe->book1D("jet_pt",   "Jet with highest p_{T} (from "+theJetCollectionLabel_.label()+");p_{T}(1^{st} jet) (GeV)",    20, 0., thePlotPhotonMaxEt_);
00150   h_jet_eta             = theDbe->book1D("jet_eta", "Jet with highest p_{T} (from "+theJetCollectionLabel_.label()+");#eta(1^{st} jet)", 20, -thePlotJetMaxEta_, thePlotJetMaxEta_);
00151   h_deltaPhi_photon_jet = theDbe->book1D("deltaPhi_photon_jet", "#Delta#phi between Highest E_{T} #gamma and jet;#Delta#phi(#gamma,1^{st} jet)", 20, 0, 3.1415926);
00152   h_deltaPhi_jet_jet2   = theDbe->book1D("deltaPhi_jet_jet2", "#Delta#phi between Highest E_{T} jet and 2^{nd} jet;#Delta#phi(1^{st} jet,2^{nd} jet)", 20, 0, 3.1415926);
00153   h_deltaEt_photon_jet  = theDbe->book1D("deltaEt_photon_jet",  "(E_{T}(#gamma)-p_{T}(jet))/E_{T}(#gamma) when #Delta#phi(#gamma,1^{st} jet) > 2.8;#DeltaE_{T}(#gamma,1^{st} jet)/E_{T}(#gamma)", 20, -1.0, 1.0);
00154   h_jet_count           = theDbe->book1D("jet_count", "Number of "+theJetCollectionLabel_.label()+" (p_{T} > "+aString+" GeV);Number of Jets", 8, -0.5, 7.5);
00155   h_jet2_pt             = theDbe->book1D("jet2_pt",   "Jet with 2^{nd} highest p_{T} (from "+theJetCollectionLabel_.label()+");p_{T}(2^{nd} jet) (GeV)",    20, 0., thePlotPhotonMaxEt_);
00156   h_jet2_eta            = theDbe->book1D("jet2_eta", "Jet with 2^{nd} highest p_{T} (from "+theJetCollectionLabel_.label()+");#eta(2^{nd} jet)", 20, -thePlotJetMaxEta_, thePlotJetMaxEta_);
00157   h_jet2_ptOverPhotonEt = theDbe->book1D("jet2_ptOverPhotonEt", "p_{T}(2^{nd} highest jet) / E_{T}(#gamma);p_{T}(2^{nd} Jet)/E_{T}(#gamma)", 20, 0.0, 4.0);
00158   h_deltaPhi_photon_jet2 = theDbe->book1D("deltaPhi_photon_jet2","#Delta#phi between Highest E_{T} #gamma and 2^{nd} highest jet;#Delta#phi(#gamma,2^{nd} jet)", 20, 0, 3.1415926);
00159   h_deltaR_jet_jet2      = theDbe->book1D("deltaR_jet_jet2", "#DeltaR between Highest Jet and 2^{nd} Highest;#DeltaR(1^{st} jet,2^{nd} jet)", 30, 0, 6.0);
00160   h_deltaR_photon_jet2   = theDbe->book1D("deltaR_photon_jet2", "#DeltaR between Highest E_{T} #gamma and 2^{nd} jet;#DeltaR(#gamma, 2^{nd} jet)", 30, 0, 6.0);
00161 
00162   // Photon Et for different jet configurations
00163   Float_t bins_et[] = {15,20,30,50,80};
00164   int num_bins_et = 4;
00165   h_photon_et_jetcs = theDbe->book1D("photon_et_jetcs", "#gamma with highest E_{T} (#eta(jet)<1.45, #eta(#gamma)#eta(jet)>0);E_{T}(#gamma) (GeV)", num_bins_et, bins_et);
00166   h_photon_et_jetco = theDbe->book1D("photon_et_jetco", "#gamma with highest E_{T} (#eta(jet)<1.45, #eta(#gamma)#eta(jet)<0);E_{T}(#gamma) (GeV)", num_bins_et, bins_et);
00167   h_photon_et_jetfs = theDbe->book1D("photon_et_jetfs", "#gamma with highest E_{T} (1.55<#eta(jet)<2.5, #eta(#gamma)#eta(jet)>0);E_{T}(#gamma) (GeV)", num_bins_et, bins_et);
00168   h_photon_et_jetfo = theDbe->book1D("photon_et_jetfo", "#gamma with highest E_{T} (1.55<#eta(jet)<2.5, #eta(#gamma)#eta(jet)<0);E_{T}(#gamma) (GeV)", num_bins_et, bins_et);
00169   h_photon_et_jetcs->getTH1F()->Sumw2();
00170   h_photon_et_jetco->getTH1F()->Sumw2();
00171   h_photon_et_jetfs->getTH1F()->Sumw2();
00172   h_photon_et_jetfo->getTH1F()->Sumw2();
00173   // Ratio of the above Photon Et distributions
00174   h_photon_et_ratio_co_cs = theDbe->book1D("photon_et_ratio_00_co_cs", "D(|#eta(jet)|<1.45, #eta(jet)*#eta(#gamma)<0) / D(|#eta(jet)|<1.45, #eta(jet)*#eta(#gamma)>0);E_{T}(#gamma) (GeV); ratio", num_bins_et, bins_et);
00175   h_photon_et_ratio_fo_fs = theDbe->book1D("photon_et_ratio_01_fo_fs", "D(1.55<|#eta(jet)|<2.6, #eta(jet)*#eta(#gamma)<0) / D(1.55<|#eta(jet)|<2.6, #eta(jet)*#eta(#gamma)>0);E_{T}(#gamma) (GeV); ratio", num_bins_et, bins_et);
00176   h_photon_et_ratio_cs_fs = theDbe->book1D("photon_et_ratio_02_cs_fs", "D(|#eta(jet)|<1.45, #eta(jet)*#eta(#gamma)>0) / D(1.55<|#eta(jet)|<2.6, #eta(jet)*#eta(#gamma)>0);E_{T}(#gamma) (GeV); ratio", num_bins_et, bins_et);
00177   h_photon_et_ratio_co_fs = theDbe->book1D("photon_et_ratio_03_co_fs", "D(|#eta(jet)|<1.45, #eta(jet)*#eta(#gamma)<0) / D(1.55<|#eta(jet)|<2.6, #eta(jet)*#eta(#gamma)>0);E_{T}(#gamma) (GeV); ratio", num_bins_et, bins_et);
00178   h_photon_et_ratio_cs_fo = theDbe->book1D("photon_et_ratio_04_cs_fo", "D(|#eta(jet)|<1.45, #eta(jet)*#eta(#gamma)>0) / D(1.55<|#eta(jet)|<2.6, #eta(jet)*#eta(#gamma)<0);E_{T}(#gamma) (GeV); ratio", num_bins_et, bins_et);
00179   h_photon_et_ratio_co_fo = theDbe->book1D("photon_et_ratio_05_co_fo", "D(|#eta(jet)|<1.45, #eta(jet)*#eta(#gamma)<0) / D(1.55<|#eta(jet)|<2.6, #eta(jet)*#eta(#gamma)<0);E_{T}(#gamma) (GeV); ratio", num_bins_et, bins_et);
00180   h_photon_et_ratio_co_cs->getTH1F()->Sumw2();
00181   h_photon_et_ratio_fo_fs->getTH1F()->Sumw2();
00182   h_photon_et_ratio_cs_fs->getTH1F()->Sumw2();
00183   h_photon_et_ratio_co_fs->getTH1F()->Sumw2();
00184   h_photon_et_ratio_cs_fo->getTH1F()->Sumw2();
00185   h_photon_et_ratio_co_fo->getTH1F()->Sumw2();
00186 }
00187 
00188 
00192 void QcdPhotonsDQM::beginRun( const edm::Run &iRun, const edm::EventSetup &iSet ) {
00193 
00194   // passed as parameter to HLTConfigProvider::init(), not yet used
00195   bool isConfigChanged = false;
00196 
00197   // isValidHltConfig_ could be used to short-circuit analyze() in case of problems
00198   isValidHltConfig_ = hltConfigProvider_.init( iRun, iSet, "HLT", isConfigChanged );
00199 
00200   num_events_in_run = 0;
00201 }
00202 
00203 
00204 void QcdPhotonsDQM::analyze(const Event& iEvent, const EventSetup& iSetup) {
00205   num_events_in_run++;
00206 
00208   //if( ! isValidHltConfig_ ) return;
00209   
00210   LogTrace(logTraceName)<<"Analysis of event # ";
00211 
00213   // Did event pass HLT paths?
00214   Handle<TriggerResults> HLTresults;
00215   iEvent.getByLabel(trigTag_, HLTresults); 
00216   if (!HLTresults.isValid()) {
00217     //LogWarning("") << ">>> TRIGGER collection does not exist !!!";
00218     return;
00219   }
00220   const edm::TriggerNames & trigNames = iEvent.triggerNames(*HLTresults);
00221 
00222   bool passed_HLT=false;
00223 
00224   // See if event passed trigger paths
00225   //  increment that bin in the trigger plot
00226   for (unsigned int i=0; i<thePlotTheseTriggersToo_.size(); i++) {
00227     passed_HLT = false;
00228     for (unsigned int ti=0; (ti<trigNames.size()) && !passed_HLT; ++ti) {
00229       size_t pos = trigNames.triggerName(ti).find(thePlotTheseTriggersToo_[i]);
00230       if (pos==0) passed_HLT = HLTresults->accept(ti);
00231     }
00232     if (passed_HLT) h_triggers_passed->Fill(i);
00233   }
00234 
00235 
00236    // grab photons
00237   Handle<PhotonCollection> photonCollection;
00238   iEvent.getByLabel(thePhotonCollectionLabel_, photonCollection);
00239 
00240   // If photon collection is empty, exit
00241   if (!photonCollection.isValid()) return;
00242 
00243 
00244   // Quit if the event did not pass the HLT path we care about
00245   passed_HLT = false;
00246   {
00247     //bool found=false;
00248   for (unsigned int ti=0; ti<trigNames.size(); ++ti) {
00249     size_t pos = trigNames.triggerName(ti).find(theTriggerPathToPass_);
00250     if (pos==0) {
00251       passed_HLT = HLTresults->accept(ti);
00252       //found=true;
00253       break;
00254     }
00255   }
00256 
00257   // Assumption: reco photons are ordered by Et
00258   for (PhotonCollection::const_iterator recoPhoton = photonCollection->begin(); recoPhoton!=photonCollection->end(); recoPhoton++){
00259     // stop looping over photons once we get to too low Et
00260     if ( recoPhoton->et() < theMinPhotonEt_ ) break;
00261 
00262     h_photon_et_beforeCuts->Fill(recoPhoton->et());
00263     break; // leading photon only
00264   }
00265 
00266   if (!passed_HLT) {
00267     return;
00268   }
00269   }
00270   
00272 
00273   //std::cout << "\tpassed main trigger (" << theTriggerPathToPass_ << ")" << std::endl;
00274 
00276   // Does event have valid vertex?
00277   // Get the primary event vertex
00278   Handle<VertexCollection> vertexHandle;
00279   iEvent.getByLabel(theVertexCollectionLabel_, vertexHandle);
00280   VertexCollection vertexCollection = *(vertexHandle.product());
00281   //double vtx_ndof = -1.0;
00282   //double vtx_z    = 0.0;
00283   //bool   vtx_isFake = true;
00284   //if (vertexCollection.size()>0) {
00285   //  vtx_ndof = vertexCollection.begin()->ndof();
00286   //  vtx_z    = vertexCollection.begin()->z();
00287   //  vtx_isFake = false;
00288   //}
00289   //if (vtx_isFake || fabs(vtx_z)>15 || vtx_ndof<4) return;
00290   
00291   int nvvertex = 0;
00292   for (unsigned int i=0; i<vertexCollection.size(); ++i) {
00293     if (vertexCollection[i].isValid()) nvvertex++;
00294   }
00295   if (nvvertex==0) return;
00296 
00298 
00299   //std::cout << "\tpassed vertex selection" << std::endl;
00300 
00301 
00303   // Did the event pass certain L1 Technical Trigger bits?
00304   // It's probably beam halo
00305   //  TODO: ADD code
00307 
00308 
00309   // For finding spikes
00310   Handle<EcalRecHitCollection> EBReducedRecHits;
00311   iEvent.getByLabel(theBarrelRecHitTag, EBReducedRecHits);
00312   Handle<EcalRecHitCollection> EEReducedRecHits;
00313   iEvent.getByLabel(theEndcapRecHitTag, EEReducedRecHits); 
00314   EcalClusterLazyTools lazyTool(iEvent, iSetup, theBarrelRecHitTag, theEndcapRecHitTag);
00315 
00316 
00317   // Find the highest et "decent" photon
00318   float photon_et  = -9.0;
00319   float photon_eta = -9.0;
00320   float photon_phi = -9.0;
00321   bool  photon_passPhotonID = false;
00322   bool  found_lead_pho = false;
00323   int   photon_count_bar = 0;
00324   int   photon_count_end = 0;
00325   // Assumption: reco photons are ordered by Et
00326   for (PhotonCollection::const_iterator recoPhoton = photonCollection->begin(); recoPhoton!=photonCollection->end(); recoPhoton++){
00327 
00328     // stop looping over photons once we get to too low Et
00329     if ( recoPhoton->et() < theMinPhotonEt_ ) break;
00330 
00331     /*
00332     //  Ignore ECAL Spikes
00333     const reco::CaloClusterPtr  seed = recoPhoton->superCluster()->seed();
00334     DetId id = lazyTool.getMaximum(*seed).first; // Cluster shape variables
00335     //    float time  = -999., outOfTimeChi2 = -999., chi2 = -999.;  // UNUSED
00336     int   flags=-1, severity = -1; 
00337     const EcalRecHitCollection & rechits = ( recoPhoton->isEB() ? *EBReducedRecHits : *EEReducedRecHits); 
00338     EcalRecHitCollection::const_iterator it = rechits.find( id );
00339     if( it != rechits.end() ) {
00340       //      time = it->time(); // UNUSED
00341       //      outOfTimeChi2 = it->outOfTimeChi2(); // UNUSED
00342       //      chi2 = it->chi2(); // UNUSED
00343       flags = it->recoFlag();
00344 
00345       edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00346       iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00347       severity = sevlv->severityLevel( id, rechits);
00348     }
00349     bool isNotSpike = ((recoPhoton->isEB() && (severity!=3 && severity!=4 ) && (flags != 2) ) || recoPhoton->isEE());
00350     if (!isNotSpike) continue;  // move on to next photon
00351     // END of determining ECAL Spikes
00352     */
00353 
00354     bool pho_current_passPhotonID = false;
00355     bool pho_current_isEB = recoPhoton->isEB();
00356     bool pho_current_isEE = recoPhoton->isEE();
00357 
00358     if ( pho_current_isEB && (recoPhoton->sigmaIetaIeta() < 0.01 || recoPhoton->hadronicOverEm() < 0.05) ) {
00359       // Photon object in barrel passes photon ID
00360       pho_current_passPhotonID = true;
00361       photon_count_bar++;
00362     } else if ( pho_current_isEE && (recoPhoton->hadronicOverEm() < 0.05) ) {
00363       // Photon object in endcap passes photon ID
00364       pho_current_passPhotonID = true;
00365       photon_count_end++;
00366     }
00367 
00368     if (!found_lead_pho) {
00369       found_lead_pho = true;
00370       photon_passPhotonID = pho_current_passPhotonID;
00371       photon_et  = recoPhoton->et();
00372       photon_eta = recoPhoton->eta();
00373       photon_phi = recoPhoton->phi();
00374     }
00375   }
00376   
00377   // If user requires a photon to be found, but none is, return.
00378   //   theRequirePhotonFound should pretty much always be set to 'True'
00379   //    except when running on qcd monte carlo just to see the jets.
00380   if ( theRequirePhotonFound_ && (!photon_passPhotonID || photon_et<theMinPhotonEt_) ) return;
00381 
00382 
00384   // Find the highest et jet
00385   Handle<View<Jet> > jetCollection;
00386   iEvent.getByLabel (theJetCollectionLabel_,jetCollection);
00387   if (!jetCollection.isValid()) return;
00388 
00389   float jet_pt    = -8.0;
00390   float jet_eta   = -8.0;
00391   float jet_phi   = -8.0;
00392   int   jet_count = 0;
00393   float jet2_pt   = -9.0;
00394   float jet2_eta  = -9.0;
00395   float jet2_phi  = -9.0;
00396   // Assumption: jets are ordered by Et
00397   for (unsigned int i_jet = 0; i_jet < jetCollection->size(); i_jet++) {
00398     const Jet* jet = & jetCollection->at(i_jet);
00399 
00400     float jet_current_pt = jet->pt();
00401 
00402     // don't care about jets that overlap with the lead photon
00403     if ( deltaR(jet->eta(), jet->phi(), photon_eta, photon_phi) < 0.5 ) continue;
00404     // stop looping over jets once we get to too low Et
00405     if (jet_current_pt < theMinJetPt_) break;
00406 
00407     jet_count++;
00408     if (jet_current_pt > jet_pt) {
00409       jet2_pt  = jet_pt;  // 2nd highest jet get's et from current highest
00410       jet2_eta = jet_eta;
00411       jet2_phi = jet_phi;
00412       jet_pt   = jet_current_pt; // current highest jet gets et from the new highest
00413       jet_eta  = jet->eta();
00414       jet_phi  = jet->phi();
00415     } else if (jet_current_pt > jet2_pt) {
00416       jet2_pt  = jet_current_pt;
00417       jet2_eta = jet->eta();
00418       jet2_phi = jet->phi();
00419     }
00420   }
00422 
00423 
00425   // Fill histograms if a jet found
00426   // NOTE: if a photon was required to be found, but wasn't
00427   //        we wouldn't have made it to this point in the code
00428   if ( jet_pt > 0.0 ) {
00429 
00430     // Photon Plots
00431     h_photon_et       ->Fill( photon_et  );
00432     h_photon_eta      ->Fill( photon_eta );
00433     h_photon_count_bar->Fill( photon_count_bar );
00434     h_photon_count_end->Fill( photon_count_end );
00435 
00436     // Photon Et hists for different orientations to the jet
00437     if ( fabs(photon_eta)<1.45 && photon_passPhotonID ) {  // Lead photon is in barrel
00438       if (fabs(jet_eta)<1.45){                          //   jet is in barrel
00439         if (photon_eta*jet_eta>0) {
00440           h_photon_et_jetcs->Fill(photon_et);
00441         } else {
00442           h_photon_et_jetco->Fill(photon_et);
00443         }
00444       } else if (jet_eta>1.55 && jet_eta<2.5) {         // jet is in endcap
00445         if (photon_eta*jet_eta>0) {
00446           h_photon_et_jetfs->Fill(photon_et);
00447         } else {
00448           h_photon_et_jetfo->Fill(photon_et);
00449         }
00450       }
00451     } // END of Lead Photon is in Barrel
00452 
00453     // Jet Plots
00454     h_jet_pt       ->Fill( jet_pt     );
00455     h_jet_eta      ->Fill( jet_eta    );
00456     h_jet_count    ->Fill( jet_count  );
00457     h_deltaPhi_photon_jet   ->Fill( abs(deltaPhi(photon_phi, jet_phi)) );
00458     if ( abs(deltaPhi(photon_phi,jet_phi))>2.8 ) h_deltaEt_photon_jet->Fill( (photon_et-jet_pt)/photon_et );
00459 
00460     // 2nd Highest Jet Plots
00461     if ( jet2_pt  > 0.0 ) {
00462       h_jet2_pt             ->Fill( jet2_pt  );
00463       h_jet2_eta            ->Fill( jet2_eta );
00464       h_jet2_ptOverPhotonEt ->Fill( jet2_pt/photon_et );
00465       h_deltaPhi_photon_jet2->Fill( abs(deltaPhi(photon_phi, jet2_phi)) );
00466       h_deltaPhi_jet_jet2   ->Fill( abs(deltaPhi(   jet_phi, jet2_phi)) );
00467       h_deltaR_jet_jet2     ->Fill( deltaR(   jet_eta,    jet_phi, jet2_eta, jet2_phi) );
00468       h_deltaR_photon_jet2  ->Fill( deltaR(photon_eta, photon_phi, jet2_eta, jet2_phi) );
00469     }
00470   } 
00471   // End of Filling histograms
00473 }
00474 
00475 
00476 void QcdPhotonsDQM::endJob(void) {}
00477 
00478 void QcdPhotonsDQM::endRun(const edm::Run& run, const edm::EventSetup& es) {
00479   if (num_events_in_run>0) { 
00480     h_triggers_passed->getTH1F()->Scale(1.0/num_events_in_run);
00481   }
00482   h_photon_et_ratio_co_cs->getTH1F()->Divide( h_photon_et_jetco->getTH1F(), h_photon_et_jetcs->getTH1F() );
00483   h_photon_et_ratio_fo_fs->getTH1F()->Divide( h_photon_et_jetfo->getTH1F(), h_photon_et_jetfs->getTH1F() );
00484   h_photon_et_ratio_cs_fs->getTH1F()->Divide( h_photon_et_jetcs->getTH1F(), h_photon_et_jetfs->getTH1F() );
00485   h_photon_et_ratio_co_fs->getTH1F()->Divide( h_photon_et_jetco->getTH1F(), h_photon_et_jetfs->getTH1F() );
00486   h_photon_et_ratio_cs_fo->getTH1F()->Divide( h_photon_et_jetcs->getTH1F(), h_photon_et_jetfo->getTH1F() );
00487   h_photon_et_ratio_co_fo->getTH1F()->Divide( h_photon_et_jetco->getTH1F(), h_photon_et_jetfo->getTH1F() );
00488 }
00489
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/DQM/Physics/src/QcdPhotonsDQM.cc
