Reference Manual

Go to the documentation of this file.
00001 /*
00002  *  See header file for a description of this class.
00003  *
00004  *  $Date: 2012/04/12 15:42:57 $
00005  *  $Revision: 1.30 $
00006  *  \author F. Chlebana - Fermilab
00007  */
00008 
00009 #include "DQMOffline/JetMET/interface/JetAnalyzer.h"
00010 #include "DataFormats/Common/interface/Handle.h"
00011 
00012 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00013 #include "DataFormats/JetReco/interface/CaloJet.h"
00014 
00015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00016 
00017 
00018 #include <string>
00019 using namespace edm;
00020 
00021 // ***********************************************************
00022 JetAnalyzer::JetAnalyzer(const edm::ParameterSet& pSet) {
00023   
00024   parameters   = pSet;
00025   _leadJetFlag = 0;
00026   _JetLoPass   = 0;
00027   _JetHiPass   = 0;
00028   _ptThreshold = 20.;
00029   _asymmetryThirdJetCut = 5.;
00030   _balanceThirdJetCut   = 0.2; 
00031   _n90HitsMin =0;
00032   _fHPDMax=1.;
00033   _resEMFMin=0.;
00034   _n90HitsMinLoose =0;
00035   _fHPDMaxLoose=1.;
00036   _resEMFMinLoose=0.;
00037   _n90HitsMinTight =0;
00038   _fHPDMaxTight=1.;
00039   _resEMFMinTight=0.;
00040   _sigmaEtaMinTight=-999.;
00041   _sigmaPhiMinTight=-999.;
00042 
00043 } 
00044   
00045 // ***********************************************************
00046 JetAnalyzer::~JetAnalyzer() { }
00047 
00048 
00049 // ***********************************************************
00050 void JetAnalyzer::beginJob(DQMStore * dbe) {
00051     
00052   jetname = "jetAnalyzer";
00053   
00054   LogTrace(jetname)<<"[JetAnalyzer] Parameters initialization";
00055   dbe->setCurrentFolder("JetMET/Jet/"+_source);
00056 
00057   jetME = dbe->book1D("jetReco", "jetReco", 3, 1, 4);
00058   jetME->setBinLabel(1,"CaloJets",1);
00059 
00060   //
00061   jetID = new reco::helper::JetIDHelper(parameters.getParameter<ParameterSet>("JetIDParams"));
00062   //
00063 
00064   fillJIDPassFrac = parameters.getParameter<int>("fillJIDPassFrac");
00065   makedijetselection = parameters.getParameter<int>("makedijetselection");
00066 
00067   // monitoring of eta parameter
00068   etaBin = parameters.getParameter<int>("etaBin");
00069   etaMin = parameters.getParameter<double>("etaMin");
00070   etaMax = parameters.getParameter<double>("etaMax");
00071 
00072   // monitoring of phi paramater
00073   phiBin = parameters.getParameter<int>("phiBin");
00074   phiMin = parameters.getParameter<double>("phiMin");
00075   phiMax = parameters.getParameter<double>("phiMax");
00076 
00077   // monitoring of the transverse momentum
00078   ptBin = parameters.getParameter<int>("ptBin");
00079   ptMin = parameters.getParameter<double>("ptMin");
00080   ptMax = parameters.getParameter<double>("ptMax");
00081 
00082   // 
00083   eBin = parameters.getParameter<int>("eBin");
00084   eMin = parameters.getParameter<double>("eMin");
00085   eMax = parameters.getParameter<double>("eMax");
00086 
00087   // 
00088   pBin = parameters.getParameter<int>("pBin");
00089   pMin = parameters.getParameter<double>("pMin");
00090   pMax = parameters.getParameter<double>("pMax");
00091 
00092   //
00093   _ptThreshold = parameters.getParameter<double>("ptThreshold");
00094   _asymmetryThirdJetCut = parameters.getParameter<double>("asymmetryThirdJetCut");
00095   _balanceThirdJetCut   = parameters.getParameter<double>("balanceThirdJetCut");
00096   _n90HitsMin = parameters.getParameter<int>("n90HitsMin");
00097   _fHPDMax = parameters.getParameter<double>("fHPDMax");
00098   _resEMFMin = parameters.getParameter<double>("resEMFMin");
00099   _sigmaEtaMinTight = parameters.getParameter<double>("sigmaEtaMinTight");
00100   _sigmaPhiMinTight = parameters.getParameter<double>("sigmaPhiMinTight");
00101 
00102   _n90HitsMinLoose = parameters.getParameter<int>("n90HitsMinLoose");
00103   _fHPDMaxLoose = parameters.getParameter<double>("fHPDMaxLoose");
00104   _resEMFMinLoose = parameters.getParameter<double>("resEMFMinLoose");
00105   _n90HitsMinTight = parameters.getParameter<int>("n90HitsMinTight");
00106   _fHPDMaxTight = parameters.getParameter<double>("fHPDMaxTight");
00107   _resEMFMinTight = parameters.getParameter<double>("resEMFMinTight");
00108 
00109 
00110   // Generic jet parameters
00111   mPt           = dbe->book1D("Pt",           "pt",                 ptBin,  ptMin,  ptMax);
00112   mEta          = dbe->book1D("Eta",          "eta",               etaBin, etaMin, etaMax);
00113   mPhi          = dbe->book1D("Phi",          "phi",               phiBin, phiMin, phiMax);
00114   mConstituents = dbe->book1D("Constituents", "# of constituents",     50,      0,    100);
00115   mHFrac        = dbe->book1D("HFrac",        "HFrac",                120,   -0.1,    1.1);
00116   mEFrac        = dbe->book1D("EFrac",        "EFrac",                120,   -0.1,    1.1);
00117 
00118 
00119   // Book NPV profiles
00120   //----------------------------------------------------------------------------
00121   mPt_profile           = dbe->bookProfile("Pt_profile",           "pt",                nbinsPV, PVlow, PVup,   ptBin,  ptMin,  ptMax);
00122   mEta_profile          = dbe->bookProfile("Eta_profile",          "eta",               nbinsPV, PVlow, PVup,  etaBin, etaMin, etaMax);
00123   mPhi_profile          = dbe->bookProfile("Phi_profile",          "phi",               nbinsPV, PVlow, PVup,  phiBin, phiMin, phiMax);
00124   mConstituents_profile = dbe->bookProfile("Constituents_profile", "# of constituents", nbinsPV, PVlow, PVup,      50,      0,    100);
00125   mHFrac_profile        = dbe->bookProfile("HFrac_profile",        "HFrac",             nbinsPV, PVlow, PVup,     120,   -0.1,    1.1);
00126   mEFrac_profile        = dbe->bookProfile("EFrac_profile",        "EFrac",             nbinsPV, PVlow, PVup,     120,   -0.1,    1.1);
00127 
00128   if (makedijetselection != 1)
00129     mNJets_profile = dbe->bookProfile("NJets_profile", "number of jets", nbinsPV, PVlow, PVup, 100, 0, 100);
00130 
00131 
00132   // Set NPV profiles x-axis title
00133   //----------------------------------------------------------------------------
00134   mPt_profile          ->setAxisTitle("nvtx",1);
00135   mEta_profile         ->setAxisTitle("nvtx",1);
00136   mPhi_profile         ->setAxisTitle("nvtx",1);
00137   mConstituents_profile->setAxisTitle("nvtx",1);
00138   mHFrac_profile       ->setAxisTitle("nvtx",1);
00139   mEFrac_profile       ->setAxisTitle("nvtx",1);
00140 
00141   if (makedijetselection != 1) {
00142     mNJets_profile->setAxisTitle("nvtx",1);
00143   }
00144 
00145 
00146   //mE                       = dbe->book1D("E", "E", eBin, eMin, eMax);
00147   //mP                       = dbe->book1D("P", "P", pBin, pMin, pMax);
00148   //  mMass                    = dbe->book1D("Mass", "Mass", 100, 0, 25);
00149   //
00150   mPhiVSEta                     = dbe->book2D("PhiVSEta", "PhiVSEta", 50, etaMin, etaMax, 24, phiMin, phiMax);
00151   if(makedijetselection!=1){
00152     mPt_1                    = dbe->book1D("Pt1", "Pt1", 20, 0, 100);   
00153     mPt_2                    = dbe->book1D("Pt2", "Pt2", 60, 0, 300);   
00154     mPt_3                    = dbe->book1D("Pt3", "Pt3", 100, 0, 5000);
00155     // Low and high pt trigger paths
00156     mPt_Lo                  = dbe->book1D("Pt_Lo", "Pt (Pass Low Pt Jet Trigger)", 20, 0, 100);   
00157     //mEta_Lo                 = dbe->book1D("Eta_Lo", "Eta (Pass Low Pt Jet Trigger)", etaBin, etaMin, etaMax);
00158     mPhi_Lo                 = dbe->book1D("Phi_Lo", "Phi (Pass Low Pt Jet Trigger)", phiBin, phiMin, phiMax);
00159     
00160     mPt_Hi                  = dbe->book1D("Pt_Hi", "Pt (Pass Hi Pt Jet Trigger)", 60, 0, 300);   
00161     mEta_Hi                 = dbe->book1D("Eta_Hi", "Eta (Pass Hi Pt Jet Trigger)", etaBin, etaMin, etaMax);
00162     mPhi_Hi                 = dbe->book1D("Phi_Hi", "Phi (Pass Hi Pt Jet Trigger)", phiBin, phiMin, phiMax);
00163     mNJets                   = dbe->book1D("NJets", "number of jets", 100, 0, 100);
00164 
00165     //mPt_Barrel_Lo            = dbe->book1D("Pt_Barrel_Lo", "Pt Barrel (Pass Low Pt Jet Trigger)", 20, 0, 100);   
00166     //mPhi_Barrel_Lo           = dbe->book1D("Phi_Barrel_Lo", "Phi Barrel (Pass Low Pt Jet Trigger)", phiBin, phiMin, phiMax);
00167     mConstituents_Barrel     = dbe->book1D("Constituents_Barrel", "Constituents Barrel above", 50, 0, 100);
00168     mHFrac_Barrel            = dbe->book1D("HFrac_Barrel", "HFrac Barrel", 100, 0, 1);
00169     mEFrac_Barrel            = dbe->book1D("EFrac_Barrel", "EFrac Barrel", 110, -0.05, 1.05);
00170     
00171     //mPt_EndCap_Lo            = dbe->book1D("Pt_EndCap_Lo", "Pt EndCap (Pass Low Pt Jet Trigger)", 20, 0, 100);   
00172     //mPhi_EndCap_Lo           = dbe->book1D("Phi_EndCap_Lo", "Phi EndCap (Pass Low Pt Jet Trigger)", phiBin, phiMin, phiMax);
00173     mConstituents_EndCap     = dbe->book1D("Constituents_EndCap", "Constituents EndCap", 50, 0, 100);
00174     mHFrac_EndCap            = dbe->book1D("HFrac_Endcap", "HFrac EndCap", 100, 0, 1);
00175     mEFrac_EndCap            = dbe->book1D("EFrac_Endcap", "EFrac EndCap", 110, -0.05, 1.05);
00176     
00177     //mPt_Forward_Lo           = dbe->book1D("Pt_Forward_Lo", "Pt Forward (Pass Low Pt Jet Trigger)", 20, 0, 100);  
00178     //mPhi_Forward_Lo          = dbe->book1D("Phi_Forward_Lo", "Phi Forward (Pass Low Pt Jet Trigger)", phiBin, phiMin, phiMax);
00179     mConstituents_Forward    = dbe->book1D("Constituents_Forward", "Constituents Forward", 50, 0, 100);
00180     mHFrac_Forward           = dbe->book1D("HFrac_Forward", "HFrac Forward", 100, 0, 1);
00181     mEFrac_Forward           = dbe->book1D("EFrac_Forward", "EFrac Forward", 110, -0.05, 1.05);
00182     
00183     mPt_Barrel_Hi            = dbe->book1D("Pt_Barrel_Hi", "Pt Barrel (Pass Hi Pt Jet Trigger)", 60, 0, 300);   
00184     mPhi_Barrel_Hi           = dbe->book1D("Phi_Barrel_Hi", "Phi Barrel (Pass Hi Pt Jet Trigger)", phiBin, phiMin, phiMax);
00185     //mConstituents_Barrel_Hi  = dbe->book1D("Constituents_Barrel_Hi", "Constituents Barrel (Pass Hi Pt Jet Trigger)", 50, 0, 100);
00186     //mHFrac_Barrel_Hi         = dbe->book1D("HFrac_Barrel_Hi", "HFrac Barrel (Pass Hi Pt Jet Trigger)", 100, 0, 1);
00187     
00188     mPt_EndCap_Hi            = dbe->book1D("Pt_EndCap_Hi", "Pt EndCap (Pass Hi Pt Jet Trigger)", 60, 0, 300);  
00189     mPhi_EndCap_Hi           = dbe->book1D("Phi_EndCap_Hi", "Phi EndCap (Pass Hi Pt Jet Trigger)", phiBin, phiMin, phiMax);
00190     //mConstituents_EndCap_Hi  = dbe->book1D("Constituents_EndCap_Hi", "Constituents EndCap (Pass Hi Pt Jet Trigger)", 50, 0, 100);
00191     //mHFrac_EndCap_Hi         = dbe->book1D("HFrac_EndCap_Hi", "HFrac EndCap (Pass Hi Pt Jet Trigger)", 100, 0, 1);
00192     
00193     mPt_Forward_Hi           = dbe->book1D("Pt_Forward_Hi", "Pt Forward (Pass Hi Pt Jet Trigger)", 60, 0, 300);  
00194     mPhi_Forward_Hi          = dbe->book1D("Phi_Forward_Hi", "Phi Forward (Pass Hi Pt Jet Trigger)", phiBin, phiMin, phiMax);
00195     //mConstituents_Forward_Hi = dbe->book1D("Constituents_Forward_Hi", "Constituents Forward (Pass Hi Pt Jet Trigger)", 50, 0, 100);
00196     //mHFrac_Forward_Hi        = dbe->book1D("HFrac_Forward_Hi", "HFrac Forward (Pass Hi Pt Jet Trigger)", 100, 0, 1);
00197     
00198     mPhi_Barrel              = dbe->book1D("Phi_Barrel", "Phi_Barrel", phiBin, phiMin, phiMax);
00199     //mE_Barrel                = dbe->book1D("E_Barrel", "E_Barrel", eBin, eMin, eMax);
00200     mPt_Barrel               = dbe->book1D("Pt_Barrel", "Pt_Barrel", ptBin, ptMin, ptMax);
00201     
00202     mPhi_EndCap              = dbe->book1D("Phi_EndCap", "Phi_EndCap", phiBin, phiMin, phiMax);
00203     //mE_EndCap                = dbe->book1D("E_EndCap", "E_EndCap", eBin, eMin, 2*eMax);
00204     mPt_EndCap               = dbe->book1D("Pt_EndCap", "Pt_EndCap", ptBin, ptMin, ptMax);
00205     
00206     mPhi_Forward             = dbe->book1D("Phi_Forward", "Phi_Forward", phiBin, phiMin, phiMax);
00207     //mE_Forward               = dbe->book1D("E_Forward", "E_Forward", eBin, eMin, 4*eMax);
00208     mPt_Forward              = dbe->book1D("Pt_Forward", "Pt_Forward", ptBin, ptMin, ptMax);
00209     
00210     // Leading Jet Parameters
00211     mEtaFirst                = dbe->book1D("EtaFirst", "EtaFirst", 100, -5, 5);
00212     mPhiFirst                = dbe->book1D("PhiFirst", "PhiFirst", 70, -3.5, 3.5);
00213     //mEFirst                  = dbe->book1D("EFirst", "EFirst", 100, 0, 1000);
00214     mPtFirst                 = dbe->book1D("PtFirst", "PtFirst", 100, 0, 500);
00215     if(fillJIDPassFrac==1){//fillJIDPassFrac defines a collection of cleaned jets, for which we will want to fill the cleaning passing fraction
00216       mLooseJIDPassFractionVSeta      = dbe->bookProfile("LooseJIDPassFractionVSeta","LooseJIDPassFractionVSeta",etaBin, etaMin, etaMax,0.,1.2);
00217       mLooseJIDPassFractionVSpt       = dbe->bookProfile("LooseJIDPassFractionVSpt","LooseJIDPassFractionVSpt",ptBin, ptMin, ptMax,0.,1.2);
00218       mTightJIDPassFractionVSeta      = dbe->bookProfile("TightJIDPassFractionVSeta","TightJIDPassFractionVSeta",etaBin, etaMin, etaMax,0.,1.2);
00219       mTightJIDPassFractionVSpt       = dbe->bookProfile("TightJIDPassFractionVSpt","TightJIDPassFractionVSpt",ptBin, ptMin, ptMax,0.,1.2);
00220 
00221 
00222     }
00223   }
00224   // CaloJet specific
00225   mMaxEInEmTowers         = dbe->book1D("MaxEInEmTowers", "MaxEInEmTowers", 100, 0, 100);
00226   mMaxEInHadTowers        = dbe->book1D("MaxEInHadTowers", "MaxEInHadTowers", 100, 0, 100);
00227   if(makedijetselection!=1) {
00228     mHadEnergyInHO          = dbe->book1D("HadEnergyInHO", "HadEnergyInHO", 100, 0, 10);
00229     mHadEnergyInHB          = dbe->book1D("HadEnergyInHB", "HadEnergyInHB", 100, 0, 50);
00230     mHadEnergyInHF          = dbe->book1D("HadEnergyInHF", "HadEnergyInHF", 100, 0, 50);
00231     mHadEnergyInHE          = dbe->book1D("HadEnergyInHE", "HadEnergyInHE", 100, 0, 100);
00232     mEmEnergyInEB           = dbe->book1D("EmEnergyInEB", "EmEnergyInEB", 100, 0, 50);
00233     mEmEnergyInEE           = dbe->book1D("EmEnergyInEE", "EmEnergyInEE", 100, 0, 50);
00234     mEmEnergyInHF           = dbe->book1D("EmEnergyInHF", "EmEnergyInHF", 120, -20, 100);
00235   }
00236   mDPhi                   = dbe->book1D("DPhi", "dPhi btw the two leading jets", 100, 0., acos(-1.));
00237   
00238   //JetID variables
00239   
00240   mresEMF                 = dbe->book1D("resEMF", "resEMF", 50, 0., 1.);
00241   mN90Hits                = dbe->book1D("N90Hits", "N90Hits", 50, 0., 50);
00242   mfHPD                   = dbe->book1D("fHPD", "fHPD", 50, 0., 1.);
00243   mfRBX                   = dbe->book1D("fRBX", "fRBX", 50, 0., 1.);
00244 
00245   //  msigmaEta                   = dbe->book1D("sigmaEta", "sigmaEta", 50, 0., 1.);
00246   //  msigmaPhi                   = dbe->book1D("sigmaPhi", "sigmaPhi", 50, 0., 0.5);
00247   
00248   if(makedijetselection==1) {
00249     mDijetAsymmetry                   = dbe->book1D("DijetAsymmetry", "DijetAsymmetry", 100, -1., 1.);
00250     mDijetBalance                     = dbe->book1D("DijetBalance",   "DijetBalance",   100, -2., 2.);
00251     if (fillJIDPassFrac==1) {
00252       mLooseJIDPassFractionVSeta  = dbe->bookProfile("LooseJIDPassFractionVSeta","LooseJIDPassFractionVSeta",50, -3., 3.,0.,1.2);
00253       mLooseJIDPassFractionVSpt   = dbe->bookProfile("LooseJIDPassFractionVSpt","LooseJIDPassFractionVSpt",ptBin, ptMin, ptMax,0.,1.2);
00254       mTightJIDPassFractionVSeta  = dbe->bookProfile("TightJIDPassFractionVSeta","TightJIDPassFractionVSeta",50, -3., 3.,0.,1.2);
00255       mTightJIDPassFractionVSpt   = dbe->bookProfile("TightJIDPassFractionVSpt","TightJIDPassFractionVSpt",ptBin, ptMin, ptMax,0.,1.2);
00256     }
00257   }
00258 }
00259 
00260 
00261 //void JetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
00262 //                        const edm::TriggerResults& triggerResults,
00263 //                        const reco::CaloJet& jet) {
00264 
00265 
00266 // ***********************************************************
00267 void JetAnalyzer::endJob() {
00268   delete jetID;
00269 }
00270 
00271 
00272 // ***********************************************************
00273 void JetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, 
00274                           const reco::CaloJetCollection& caloJets,
00275                           const int numPV) {
00276   int numofjets=0;
00277   double  fstPhi=0.;
00278   double  sndPhi=0.;
00279   double  diff = 0.;
00280   double  corr = 0.;
00281   double  dphi = -999. ;
00282   bool thiscleaned=false;
00283   bool Loosecleaned=false;
00284   bool Tightcleaned=false;
00285   bool thisemfclean=true;
00286   bool emfcleanLoose=true;
00287   bool emfcleanTight=true;
00288 
00289   srand( iEvent.id().event() % 10000);
00290 
00291 
00292   if (makedijetselection == 1) {
00293     //Dijet selection - careful: the pT is uncorrected!
00294     //if(makedijetselection==1 && caloJets.size()>=2){
00295     if(caloJets.size()>=2){
00296       double  dphiDJ = -999. ;
00297       bool emfcleanLooseFirstJet=true;
00298       bool emfcleanLooseSecondJet=true;
00299       bool emfcleanTightFirstJet=true;
00300       bool emfcleanTightSecondJet=true;
00301       bool LoosecleanedFirstJet = false;
00302       bool LoosecleanedSecondJet = false;
00303       bool TightcleanedFirstJet = false;
00304       bool TightcleanedSecondJet = false;
00305       //both jets pass pt threshold
00306       if ((caloJets.at(0)).pt() > _ptThreshold && (caloJets.at(1)).pt() > _ptThreshold ) {
00307         if(fabs((caloJets.at(0)).eta())<3. && fabs((caloJets.at(1)).eta())<3. ){
00308           //calculate dphi
00309           dphiDJ = fabs((caloJets.at(0)).phi()-(caloJets.at(1)).phi());
00310           if (dphiDJ > 3.14) dphiDJ=fabs(dphiDJ -6.28 );
00311           //fill DPhi histo (before cutting)
00312           if (mDPhi) mDPhi->Fill (dphiDJ);
00313           //dphi cut
00314           if(fabs(dphiDJ)>2.1){
00315             //JetID 
00316             emfcleanLooseFirstJet=true;
00317             emfcleanTightFirstJet=true;
00318             emfcleanLooseSecondJet=true;
00319             emfcleanTightSecondJet=true;
00320             //jetID for first jet
00321             jetID->calculate(iEvent, (caloJets.at(0)));
00322             if(jetID->restrictedEMF()<_resEMFMinLoose && fabs((caloJets.at(0)).eta())<2.6) emfcleanLooseFirstJet=false;
00323             if(jetID->restrictedEMF()<_resEMFMinTight && fabs((caloJets.at(0)).eta())<2.6) emfcleanTightFirstJet=false;
00324             if(jetID->n90Hits()>=_n90HitsMinLoose && jetID->fHPD()<_fHPDMaxLoose && emfcleanLooseFirstJet) LoosecleanedFirstJet=true;
00325             if(jetID->n90Hits()>=_n90HitsMinTight && jetID->fHPD()<_fHPDMaxTight && sqrt((caloJets.at(0)).etaetaMoment())>_sigmaEtaMinTight && sqrt((caloJets.at(0)).phiphiMoment())>_sigmaPhiMinTight && emfcleanTightFirstJet) TightcleanedFirstJet=true;
00326             //fill the JID variables histograms BEFORE you cut on them
00327             if (mN90Hits)         mN90Hits->Fill (jetID->n90Hits());
00328             if (mfHPD)            mfHPD->Fill (jetID->fHPD());
00329             if (mresEMF)         mresEMF->Fill (jetID->restrictedEMF());
00330             if (mfRBX)            mfRBX->Fill (jetID->fRBX());
00331 
00332             //jetID for second jet
00333             jetID->calculate(iEvent, (caloJets.at(1)));
00334             if(jetID->restrictedEMF()<_resEMFMinLoose && fabs((caloJets.at(1)).eta())<2.6) emfcleanLooseSecondJet=false;
00335             if(jetID->restrictedEMF()<_resEMFMinTight && fabs((caloJets.at(1)).eta())<2.6) emfcleanTightSecondJet=false;
00336             if(jetID->n90Hits()>=_n90HitsMinLoose && jetID->fHPD()<_fHPDMaxLoose && emfcleanLooseSecondJet) LoosecleanedSecondJet=true;
00337             if(jetID->n90Hits()>=_n90HitsMinTight && jetID->fHPD()<_fHPDMaxTight && sqrt((caloJets.at(1)).etaetaMoment())>_sigmaEtaMinTight && sqrt((caloJets.at(1)).phiphiMoment())>_sigmaPhiMinTight && emfcleanTightSecondJet) TightcleanedSecondJet=true;
00338             //fill the JID variables histograms BEFORE you cut on them
00339             if (mN90Hits)         mN90Hits->Fill (jetID->n90Hits());
00340             if (mfHPD)            mfHPD->Fill (jetID->fHPD());
00341             if (mresEMF)         mresEMF->Fill (jetID->restrictedEMF());
00342             if (mfRBX)            mfRBX->Fill (jetID->fRBX());
00343 
00344             if(LoosecleanedFirstJet && LoosecleanedSecondJet) { //only if both jets are (loose) cleaned
00345               //fill histos for first jet
00346               if (mPt)   mPt->Fill ((caloJets.at(0)).pt());
00347               if (mEta)  mEta->Fill ((caloJets.at(0)).eta());
00348               if (mPhi)  mPhi->Fill ((caloJets.at(0)).phi());
00349               if (mPhiVSEta) mPhiVSEta->Fill((caloJets.at(0)).eta(),(caloJets.at(0)).phi());
00350               if (mConstituents) mConstituents->Fill ((caloJets.at(0)).nConstituents());
00351               if (mHFrac)        mHFrac->Fill ((caloJets.at(0)).energyFractionHadronic());
00352               if (mEFrac)        mEFrac->Fill ((caloJets.at(0)).emEnergyFraction());
00353               //if (mE)    mE->Fill ((caloJets.at(0)).energy());
00354               //if (mP)    mP->Fill ((caloJets.at(0)).p());
00355               //if (mMass) mMass->Fill ((caloJets.at(0)).mass());
00356               if (mMaxEInEmTowers)  mMaxEInEmTowers->Fill ((caloJets.at(0)).maxEInEmTowers());
00357               if (mMaxEInHadTowers) mMaxEInHadTowers->Fill ((caloJets.at(0)).maxEInHadTowers());
00358               if (mN90Hits)         mN90Hits->Fill (jetID->n90Hits());
00359               if (mfHPD)            mfHPD->Fill (jetID->fHPD());
00360               if (mresEMF)         mresEMF->Fill (jetID->restrictedEMF());
00361               if (mfRBX)            mfRBX->Fill (jetID->fRBX());
00362               //sigmaeta and sigmaphi only used in the tight selection.
00363               //fill the histos for them AFTER the loose selection 
00364               //  if (msigmaEta)  msigmaEta->Fill(sqrt((caloJets.at(0)).etaetaMoment()));
00365               //  if (msigmaPhi)  msigmaPhi->Fill(sqrt((caloJets.at(0)).phiphiMoment()));
00366               //fill histos for second jet
00367               if (mPt)   mPt->Fill ((caloJets.at(1)).pt());
00368               if (mEta)  mEta->Fill ((caloJets.at(1)).eta());
00369               if (mPhi)  mPhi->Fill ((caloJets.at(1)).phi());
00370               if (mPhiVSEta) mPhiVSEta->Fill((caloJets.at(1)).eta(),(caloJets.at(1)).phi());
00371               if (mConstituents) mConstituents->Fill ((caloJets.at(1)).nConstituents());
00372               if (mHFrac)        mHFrac->Fill ((caloJets.at(1)).energyFractionHadronic());
00373               if (mEFrac)        mEFrac->Fill ((caloJets.at(1)).emEnergyFraction());
00374               //if (mE)    mE->Fill ((caloJets.at(1)).energy());
00375               //if (mP)    mP->Fill ((caloJets.at(1)).p());
00376               //if (mMass) mMass->Fill ((caloJets.at(1)).mass());
00377               if (mMaxEInEmTowers)  mMaxEInEmTowers->Fill ((caloJets.at(1)).maxEInEmTowers());
00378               if (mMaxEInHadTowers) mMaxEInHadTowers->Fill ((caloJets.at(1)).maxEInHadTowers());
00379               //sigmaeta and sigmaphi only used in the tight selection.
00380               //fill the histos for them AFTER the loose selection 
00381               //  if (msigmaEta)  msigmaEta->Fill(sqrt((caloJets.at(1)).etaetaMoment()));
00382               //  if (msigmaPhi)  msigmaPhi->Fill(sqrt((caloJets.at(1)).phiphiMoment()));
00383 
00384 
00385               // Fill NPV profiles
00386               //----------------------------------------------------------------
00387               for (int ijet=0; ijet<2; ijet++) {
00388 
00389                 if (mPt_profile)           mPt_profile          ->Fill(numPV, (caloJets.at(ijet)).pt());
00390                 if (mEta_profile)          mEta_profile         ->Fill(numPV, (caloJets.at(ijet)).eta());
00391                 if (mPhi_profile)          mPhi_profile         ->Fill(numPV, (caloJets.at(ijet)).phi());
00392                 if (mConstituents_profile) mConstituents_profile->Fill(numPV, (caloJets.at(ijet)).nConstituents());
00393                 if (mHFrac_profile)        mHFrac_profile       ->Fill(numPV, (caloJets.at(ijet)).energyFractionHadronic());
00394                 if (mEFrac_profile)        mEFrac_profile       ->Fill(numPV, (caloJets.at(ijet)).emEnergyFraction());
00395               }
00396             }
00397 
00398 
00399             //let's see how many of these jets passed the JetID cleaning
00400             if(fillJIDPassFrac==1) {
00401               if(LoosecleanedFirstJet) {
00402                 mLooseJIDPassFractionVSeta->Fill((caloJets.at(0)).eta(),1.);
00403                 mLooseJIDPassFractionVSpt->Fill((caloJets.at(0)).pt(),1.);
00404               } else  {
00405                 mLooseJIDPassFractionVSeta->Fill((caloJets.at(0)).eta(),0.);
00406                 mLooseJIDPassFractionVSpt->Fill((caloJets.at(0)).pt(),0.);
00407               }
00408               if(LoosecleanedSecondJet) {
00409                 mLooseJIDPassFractionVSeta->Fill((caloJets.at(1)).eta(),1.);
00410                 mLooseJIDPassFractionVSpt->Fill((caloJets.at(1)).pt(),1.);
00411               } else  {
00412                 mLooseJIDPassFractionVSeta->Fill((caloJets.at(1)).eta(),0.);
00413                 mLooseJIDPassFractionVSpt->Fill((caloJets.at(1)).pt(),0.);
00414               }
00415               //TIGHT JID
00416               if(TightcleanedFirstJet) {
00417                 mTightJIDPassFractionVSeta->Fill((caloJets.at(0)).eta(),1.);
00418                 mTightJIDPassFractionVSpt->Fill((caloJets.at(0)).pt(),1.);
00419               } else  {
00420                 mTightJIDPassFractionVSeta->Fill((caloJets.at(0)).eta(),0.);
00421                 mTightJIDPassFractionVSpt->Fill((caloJets.at(0)).pt(),0.);
00422               }
00423               if(TightcleanedSecondJet) {
00424                 mTightJIDPassFractionVSeta->Fill((caloJets.at(1)).eta(),1.);
00425                 mTightJIDPassFractionVSpt->Fill((caloJets.at(1)).pt(),1.);
00426               } else  {
00427                 mTightJIDPassFractionVSeta->Fill((caloJets.at(1)).eta(),0.);
00428                 mTightJIDPassFractionVSpt->Fill((caloJets.at(1)).pt(),0.);
00429               }
00430 
00431             }//if fillJIDPassFrac
00432           }// FABS DPHI < 2.1
00433         }// fabs eta < 3
00434       }// pt jets > threshold
00435       //now do the dijet balance and asymmetry calculations
00436       if (fabs(caloJets.at(0).eta() < 1.4)) {
00437         double pt_dijet = (caloJets.at(0).pt() + caloJets.at(1).pt())/2;
00438         
00439         double dPhi = fabs((caloJets.at(0)).phi()-(caloJets.at(1)).phi());
00440         if (dPhi > 3.14) dPhi=fabs(dPhi -6.28 );
00441         
00442         if (dPhi > 2.7) {
00443           double pt_probe;
00444           double pt_barrel;
00445           int jet1, jet2;
00446 
00447           int randJet = rand() % 2;
00448 
00449           if (fabs(caloJets.at(1).eta() < 1.4)) {
00450             if (randJet) {
00451               jet1 = 0;
00452               jet2 = 1;
00453             }
00454             else {
00455               jet1 = 1;
00456               jet2 = 0;
00457             }
00458           
00459             /***Di-Jet Asymmetry****
00460              * leading jets eta < 1.4
00461              * leading jets dphi > 2.7
00462              * pt_third jet < threshold
00463              * A = (pt_1 - pt_2)/(pt_1 + pt_2)
00464              * jets 1 and two are randomly ordered
00465              */
00466             bool thirdJetCut = true;
00467             for (unsigned int third = 2; third < caloJets.size(); ++third) 
00468               if (caloJets.at(third).pt() > _asymmetryThirdJetCut) 
00469                 thirdJetCut = false;
00470             if (thirdJetCut) {
00471               double dijetAsymmetry = (caloJets.at(jet1).pt() - caloJets.at(jet2).pt()) / (caloJets.at(jet1).pt() + caloJets.at(jet2).pt());
00472               mDijetAsymmetry->Fill(dijetAsymmetry);
00473             }// end restriction on third jet pt in asymmetry calculation
00474               
00475           }
00476           else {
00477             jet1 = 0;
00478             jet2 = 1;
00479           }
00480           
00481           pt_barrel = caloJets.at(jet1).pt();
00482           pt_probe  = caloJets.at(jet2).pt();
00483           
00484           //dijet balance cuts
00485           /***Di-Jet Balance****
00486            * pt_dijet = (pt_probe+pt_barrel)/2
00487            * leading jets dphi > 2.7
00488            * reject evnets where pt_third/pt_dijet > 0.2
00489            * pv selection
00490            * B = (pt_probe - pt_barrel)/pt_dijet
00491            * select probe randomly from 2 jets if both leading jets are in the barrel
00492            */
00493           bool thirdJetCut = true;
00494           for (unsigned int third = 2; third < caloJets.size(); ++third) 
00495             if (caloJets.at(third).pt()/pt_dijet > _balanceThirdJetCut) 
00496               thirdJetCut = false;
00497           if (thirdJetCut) {
00498             double dijetBalance = (pt_probe - pt_barrel) / pt_dijet;
00499             mDijetBalance->Fill(dijetBalance);
00500           }// end restriction on third jet pt ratio in balance calculation
00501         }// dPhi > 2.7
00502       }// leading jet eta cut for asymmetry and balance calculations
00503     }//jet size >= 2
00504   }// do dijet selection
00505   else {
00506     for (reco::CaloJetCollection::const_iterator jet = caloJets.begin(); jet!=caloJets.end(); ++jet) {
00507       LogTrace(jetname)<<"[JetAnalyzer] Analyze Calo Jet";
00508       Loosecleaned=false;
00509       Tightcleaned=false;
00510       if (jet == caloJets.begin()) {
00511         fstPhi = jet->phi();
00512         _leadJetFlag = 1;
00513       } else {
00514         _leadJetFlag = 0;
00515       }
00516       if (jet == (caloJets.begin()+1)) sndPhi = jet->phi();
00517       //jetID
00518       jetID->calculate(iEvent, *jet);
00519       //minimal (uncorrected!) pT cut
00520       if (jet->pt() > _ptThreshold) {
00521         //  if (msigmaEta)  msigmaEta->Fill(sqrt(jet->etaetaMoment()));
00522         //  if (msigmaPhi)  msigmaPhi->Fill(sqrt(jet->phiphiMoment()));
00523         //cleaning to use for filling histograms
00524         thisemfclean=true;
00525         if(jetID->restrictedEMF()<_resEMFMin && fabs(jet->eta())<2.6) thisemfclean=false;
00526         if(jetID->n90Hits()>=_n90HitsMin && jetID->fHPD()<_fHPDMax && thisemfclean) thiscleaned=true;
00527         //loose and tight cleaning, used to fill the JetIDPAssFraction histos
00528         if(jetID->n90Hits()>=_n90HitsMinLoose && jetID->fHPD()<_fHPDMaxLoose && emfcleanLoose) Loosecleaned=true;
00529         if(jetID->n90Hits()>=_n90HitsMinTight && jetID->fHPD()<_fHPDMaxTight && sqrt(jet->etaetaMoment())>_sigmaEtaMinTight && sqrt(jet->phiphiMoment())>_sigmaPhiMinTight && emfcleanTight) Tightcleaned=true;
00530 
00531         if(fillJIDPassFrac==1) {
00532           if(Loosecleaned) {
00533             mLooseJIDPassFractionVSeta->Fill(jet->eta(),1.);
00534             mLooseJIDPassFractionVSpt->Fill(jet->pt(),1.);
00535           } else {
00536             mLooseJIDPassFractionVSeta->Fill(jet->eta(),0.);
00537             mLooseJIDPassFractionVSpt->Fill(jet->pt(),0.);
00538           }
00539           //TIGHT
00540           if(Tightcleaned) {
00541             mTightJIDPassFractionVSeta->Fill(jet->eta(),1.);
00542             mTightJIDPassFractionVSpt->Fill(jet->pt(),1.);
00543           } else {
00544             mTightJIDPassFractionVSeta->Fill(jet->eta(),0.);
00545             mTightJIDPassFractionVSpt->Fill(jet->pt(),0.);
00546           }
00547         }
00548         //eventually we could define the "cleaned" flag differently for e.g. HF
00549         if(thiscleaned) {
00550           numofjets++ ;
00551           jetME->Fill(1);      
00552           
00553           // Leading jet
00554           // Histograms are filled once per event
00555           if (_leadJetFlag == 1) { 
00556             if (mEtaFirst) mEtaFirst->Fill (jet->eta());
00557             if (mPhiFirst) mPhiFirst->Fill (jet->phi());
00558             //if (mEFirst)   mEFirst->Fill (jet->energy());
00559             if (mPtFirst)  mPtFirst->Fill (jet->pt());
00560           }
00561           // --- Passed the low pt jet trigger
00562           if (_JetLoPass == 1) {
00563           /*  if (fabs(jet->eta()) <= 1.3) {
00564               if (mPt_Barrel_Lo)           mPt_Barrel_Lo->Fill(jet->pt());
00565               if (mEta_Lo)          mEta_Lo->Fill(jet->eta());
00566               if (mPhi_Barrel_Lo)          mPhi_Barrel_Lo->Fill(jet->phi());
00567             }
00568             if ( (fabs(jet->eta()) > 1.3) && (fabs(jet->eta()) <= 3) ) {
00569               if (mPt_EndCap_Lo)           mPt_EndCap_Lo->Fill(jet->pt());
00570               if (mEta_Lo)          mEta_Lo->Fill(jet->eta());
00571               if (mPhi_EndCap_Lo)          mPhi_EndCap_Lo->Fill(jet->phi());
00572             }
00573             if (fabs(jet->eta()) > 3.0) {
00574               if (mPt_Forward_Lo)           mPt_Forward_Lo->Fill(jet->pt());
00575               if (mEta_Lo)          mEta_Lo->Fill(jet->eta());
00576               if (mPhi_Forward_Lo)          mPhi_Forward_Lo->Fill(jet->phi());
00577             } */
00578             //if (mEta_Lo) mEta_Lo->Fill (jet->eta());
00579             if (mPhi_Lo) mPhi_Lo->Fill (jet->phi());
00580             if (mPt_Lo)  mPt_Lo->Fill (jet->pt());
00581           }
00582           
00583           // --- Passed the high pt jet trigger
00584           if (_JetHiPass == 1) {
00585             if (fabs(jet->eta()) <= 1.3) {
00586               if (mPt_Barrel_Hi && jet->pt()>100.)           mPt_Barrel_Hi->Fill(jet->pt());
00587               if (mEta_Hi && jet->pt()>100.)          mEta_Hi->Fill(jet->eta());
00588               if (mPhi_Barrel_Hi)          mPhi_Barrel_Hi->Fill(jet->phi());
00589               //if (mConstituents_Barrel_Hi) mConstituents_Barrel_Hi->Fill(jet->nConstituents());       
00590               //if (mHFrac_Barrel_Hi)        mHFrac_Barrel_Hi->Fill(jet->energyFractionHadronic());     
00591             }
00592             if ( (fabs(jet->eta()) > 1.3) && (fabs(jet->eta()) <= 3) ) {
00593               if (mPt_EndCap_Hi && jet->pt()>100.)           mPt_EndCap_Hi->Fill(jet->pt());
00594               if (mEta_Hi && jet->pt()>100.)          mEta_Hi->Fill(jet->eta());
00595               if (mPhi_EndCap_Hi)          mPhi_EndCap_Hi->Fill(jet->phi());
00596               //if (mConstituents_EndCap_Hi) mConstituents_EndCap_Hi->Fill(jet->nConstituents());       
00597               //if (mHFrac_EndCap_Hi)        mHFrac_EndCap_Hi->Fill(jet->energyFractionHadronic());     
00598             }
00599             if (fabs(jet->eta()) > 3.0) {
00600               if (mPt_Forward_Hi && jet->pt()>100.)           mPt_Forward_Hi->Fill(jet->pt());
00601               if (mEta_Hi && jet->pt()>100.)          mEta_Hi->Fill(jet->eta());
00602               if (mPhi_Forward_Hi)          mPhi_Forward_Hi->Fill(jet->phi());
00603               //if (mConstituents_Forward_Hi) mConstituents_Forward_Hi->Fill(jet->nConstituents());     
00604               //if (mHFrac_Forward_Hi)        mHFrac_Forward_Hi->Fill(jet->energyFractionHadronic());   
00605             }
00606             
00607             if (mEta_Hi && jet->pt()>100.) mEta_Hi->Fill (jet->eta());
00608             if (mPhi_Hi) mPhi_Hi->Fill (jet->phi());
00609             if (mPt_Hi)  mPt_Hi->Fill (jet->pt());
00610           }
00611           
00612           if (mPt)   mPt->Fill (jet->pt());
00613           if (mPt_1) mPt_1->Fill (jet->pt());
00614           if (mPt_2) mPt_2->Fill (jet->pt());
00615           if (mPt_3) mPt_3->Fill (jet->pt());
00616           if (mEta)  mEta->Fill (jet->eta());
00617           if (mPhi)  mPhi->Fill (jet->phi());
00618           
00619           if (mPhiVSEta) mPhiVSEta->Fill(jet->eta(),jet->phi());
00620           
00621           if (mConstituents) mConstituents->Fill (jet->nConstituents());
00622           if (mHFrac)        mHFrac->Fill (jet->energyFractionHadronic());
00623           if (mEFrac)        mEFrac->Fill (jet->emEnergyFraction());
00624           
00625 
00626           // Fill NPV profiles
00627           //--------------------------------------------------------------------
00628           if (mPt_profile)           mPt_profile          ->Fill(numPV, jet->pt());
00629           if (mEta_profile)          mEta_profile         ->Fill(numPV, jet->eta());
00630           if (mPhi_profile)          mPhi_profile         ->Fill(numPV, jet->phi());
00631           if (mConstituents_profile) mConstituents_profile->Fill(numPV, jet->nConstituents());
00632           if (mHFrac_profile)        mHFrac_profile       ->Fill(numPV, jet->energyFractionHadronic());
00633           if (mEFrac_profile)        mEFrac_profile       ->Fill(numPV, jet->emEnergyFraction());
00634 
00635 
00636           if (fabs(jet->eta()) <= 1.3) {
00637             if (mPt_Barrel)   mPt_Barrel->Fill (jet->pt());
00638             if (mPhi_Barrel)  mPhi_Barrel->Fill (jet->phi());
00639             //if (mE_Barrel)    mE_Barrel->Fill (jet->energy());
00640       if (mConstituents_Barrel)    mConstituents_Barrel->Fill(jet->nConstituents());    
00641       if (mHFrac_Barrel)           mHFrac_Barrel->Fill(jet->energyFractionHadronic());  
00642       if (mEFrac_Barrel)           mEFrac_Barrel->Fill(jet->emEnergyFraction());        
00643           }
00644           if ( (fabs(jet->eta()) > 1.3) && (fabs(jet->eta()) <= 3) ) {
00645             if (mPt_EndCap)   mPt_EndCap->Fill (jet->pt());
00646             if (mPhi_EndCap)  mPhi_EndCap->Fill (jet->phi());
00647             //if (mE_EndCap)    mE_EndCap->Fill (jet->energy());
00648       if (mConstituents_EndCap)    mConstituents_EndCap->Fill(jet->nConstituents());    
00649       if (mHFrac_EndCap)           mHFrac_EndCap->Fill(jet->energyFractionHadronic());
00650       if (mEFrac_EndCap)           mEFrac_EndCap->Fill(jet->emEnergyFraction());        
00651           }
00652           if (fabs(jet->eta()) > 3.0) {
00653             if (mPt_Forward)   mPt_Forward->Fill (jet->pt());
00654             if (mPhi_Forward)  mPhi_Forward->Fill (jet->phi());
00655             //if (mE_Forward)    mE_Forward->Fill (jet->energy());
00656       if (mConstituents_Forward)    mConstituents_Forward->Fill(jet->nConstituents());  
00657       if (mHFrac_Forward)           mHFrac_Forward->Fill(jet->energyFractionHadronic());
00658       if (mEFrac_Forward)           mEFrac_Forward->Fill(jet->emEnergyFraction());      
00659           }
00660           
00661           //if (mE)    mE->Fill (jet->energy());
00662           //if (mP)    mP->Fill (jet->p());
00663           // if (mMass) mMass->Fill (jet->mass());
00664           
00665           if (mMaxEInEmTowers)  mMaxEInEmTowers->Fill (jet->maxEInEmTowers());
00666           if (mMaxEInHadTowers) mMaxEInHadTowers->Fill (jet->maxEInHadTowers());
00667           
00668           if (mHadEnergyInHO)   mHadEnergyInHO->Fill (jet->hadEnergyInHO());
00669           if (mHadEnergyInHB)   mHadEnergyInHB->Fill (jet->hadEnergyInHB());
00670           if (mHadEnergyInHF)   mHadEnergyInHF->Fill (jet->hadEnergyInHF());
00671           if (mHadEnergyInHE)   mHadEnergyInHE->Fill (jet->hadEnergyInHE());
00672           if (mEmEnergyInEB)    mEmEnergyInEB->Fill (jet->emEnergyInEB());
00673           if (mEmEnergyInEE)    mEmEnergyInEE->Fill (jet->emEnergyInEE());
00674           if (mEmEnergyInHF)    mEmEnergyInHF->Fill (jet->emEnergyInHF());
00675           
00676           if (mN90Hits)         mN90Hits->Fill (jetID->n90Hits());
00677           if (mfHPD)            mfHPD->Fill (jetID->fHPD());
00678           if (mresEMF)         mresEMF->Fill (jetID->restrictedEMF());
00679           if (mfRBX)            mfRBX->Fill (jetID->fRBX());
00680           
00681           //calculate correctly the dphi
00682           if(numofjets>1) {
00683             diff = fabs(fstPhi - sndPhi);
00684             corr = 2*acos(-1.) - diff;
00685             if(diff < acos(-1.)) { 
00686               dphi = diff; 
00687             } else { 
00688               dphi = corr;
00689             }
00690           }
00691         }
00692       }//pt cut
00693     }
00694 
00695     if (mNJets) mNJets->Fill (numofjets);
00696     if (mDPhi && dphi>-998.) mDPhi->Fill (dphi);
00697 
00698 
00699     if (mNJets_profile) mNJets_profile->Fill(numPV, numofjets);
00700 
00701 
00702   }//not dijet
00703 }
00704
/data/git/CMSSW_5_3_11_patch5/src/DQMOffline/JetMET/src/JetAnalyzer.cc
