CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
JetAnalyzer.cc
Go to the documentation of this file.
1 /*
2  * See header file for a description of this class.
3  *
4  * $Date: 2012/04/12 15:42:57 $
5  * $Revision: 1.30 $
6  * \author F. Chlebana - Fermilab
7  */
8 
11 
14 
16 
17 
18 #include <string>
19 using namespace edm;
20 
21 // ***********************************************************
23 
24  parameters = pSet;
25  _leadJetFlag = 0;
26  _JetLoPass = 0;
27  _JetHiPass = 0;
28  _ptThreshold = 20.;
29  _asymmetryThirdJetCut = 5.;
30  _balanceThirdJetCut = 0.2;
31  _n90HitsMin =0;
32  _fHPDMax=1.;
33  _resEMFMin=0.;
34  _n90HitsMinLoose =0;
35  _fHPDMaxLoose=1.;
36  _resEMFMinLoose=0.;
37  _n90HitsMinTight =0;
38  _fHPDMaxTight=1.;
39  _resEMFMinTight=0.;
40  _sigmaEtaMinTight=-999.;
41  _sigmaPhiMinTight=-999.;
42 
43 }
44 
45 // ***********************************************************
47 
48 
49 // ***********************************************************
51 
52  jetname = "jetAnalyzer";
53 
54  LogTrace(jetname)<<"[JetAnalyzer] Parameters initialization";
55  dbe->setCurrentFolder("JetMET/Jet/"+_source);
56 
57  jetME = dbe->book1D("jetReco", "jetReco", 3, 1, 4);
58  jetME->setBinLabel(1,"CaloJets",1);
59 
60  //
61  jetID = new reco::helper::JetIDHelper(parameters.getParameter<ParameterSet>("JetIDParams"));
62  //
63 
64  fillJIDPassFrac = parameters.getParameter<int>("fillJIDPassFrac");
65  makedijetselection = parameters.getParameter<int>("makedijetselection");
66 
67  // monitoring of eta parameter
68  etaBin = parameters.getParameter<int>("etaBin");
69  etaMin = parameters.getParameter<double>("etaMin");
70  etaMax = parameters.getParameter<double>("etaMax");
71 
72  // monitoring of phi paramater
73  phiBin = parameters.getParameter<int>("phiBin");
74  phiMin = parameters.getParameter<double>("phiMin");
75  phiMax = parameters.getParameter<double>("phiMax");
76 
77  // monitoring of the transverse momentum
78  ptBin = parameters.getParameter<int>("ptBin");
79  ptMin = parameters.getParameter<double>("ptMin");
80  ptMax = parameters.getParameter<double>("ptMax");
81 
82  //
83  eBin = parameters.getParameter<int>("eBin");
84  eMin = parameters.getParameter<double>("eMin");
85  eMax = parameters.getParameter<double>("eMax");
86 
87  //
88  pBin = parameters.getParameter<int>("pBin");
89  pMin = parameters.getParameter<double>("pMin");
90  pMax = parameters.getParameter<double>("pMax");
91 
92  //
93  _ptThreshold = parameters.getParameter<double>("ptThreshold");
94  _asymmetryThirdJetCut = parameters.getParameter<double>("asymmetryThirdJetCut");
95  _balanceThirdJetCut = parameters.getParameter<double>("balanceThirdJetCut");
96  _n90HitsMin = parameters.getParameter<int>("n90HitsMin");
97  _fHPDMax = parameters.getParameter<double>("fHPDMax");
98  _resEMFMin = parameters.getParameter<double>("resEMFMin");
99  _sigmaEtaMinTight = parameters.getParameter<double>("sigmaEtaMinTight");
100  _sigmaPhiMinTight = parameters.getParameter<double>("sigmaPhiMinTight");
101 
102  _n90HitsMinLoose = parameters.getParameter<int>("n90HitsMinLoose");
103  _fHPDMaxLoose = parameters.getParameter<double>("fHPDMaxLoose");
104  _resEMFMinLoose = parameters.getParameter<double>("resEMFMinLoose");
105  _n90HitsMinTight = parameters.getParameter<int>("n90HitsMinTight");
106  _fHPDMaxTight = parameters.getParameter<double>("fHPDMaxTight");
107  _resEMFMinTight = parameters.getParameter<double>("resEMFMinTight");
108 
109 
110  // Generic jet parameters
111  mPt = dbe->book1D("Pt", "pt", ptBin, ptMin, ptMax);
112  mEta = dbe->book1D("Eta", "eta", etaBin, etaMin, etaMax);
113  mPhi = dbe->book1D("Phi", "phi", phiBin, phiMin, phiMax);
114  mConstituents = dbe->book1D("Constituents", "# of constituents", 50, 0, 100);
115  mHFrac = dbe->book1D("HFrac", "HFrac", 120, -0.1, 1.1);
116  mEFrac = dbe->book1D("EFrac", "EFrac", 120, -0.1, 1.1);
117 
118 
119  // Book NPV profiles
120  //----------------------------------------------------------------------------
121  mPt_profile = dbe->bookProfile("Pt_profile", "pt", nbinsPV, PVlow, PVup, ptBin, ptMin, ptMax);
122  mEta_profile = dbe->bookProfile("Eta_profile", "eta", nbinsPV, PVlow, PVup, etaBin, etaMin, etaMax);
123  mPhi_profile = dbe->bookProfile("Phi_profile", "phi", nbinsPV, PVlow, PVup, phiBin, phiMin, phiMax);
124  mConstituents_profile = dbe->bookProfile("Constituents_profile", "# of constituents", nbinsPV, PVlow, PVup, 50, 0, 100);
125  mHFrac_profile = dbe->bookProfile("HFrac_profile", "HFrac", nbinsPV, PVlow, PVup, 120, -0.1, 1.1);
126  mEFrac_profile = dbe->bookProfile("EFrac_profile", "EFrac", nbinsPV, PVlow, PVup, 120, -0.1, 1.1);
127 
128  if (makedijetselection != 1)
129  mNJets_profile = dbe->bookProfile("NJets_profile", "number of jets", nbinsPV, PVlow, PVup, 100, 0, 100);
130 
131 
132  // Set NPV profiles x-axis title
133  //----------------------------------------------------------------------------
134  mPt_profile ->setAxisTitle("nvtx",1);
135  mEta_profile ->setAxisTitle("nvtx",1);
136  mPhi_profile ->setAxisTitle("nvtx",1);
137  mConstituents_profile->setAxisTitle("nvtx",1);
138  mHFrac_profile ->setAxisTitle("nvtx",1);
139  mEFrac_profile ->setAxisTitle("nvtx",1);
140 
141  if (makedijetselection != 1) {
142  mNJets_profile->setAxisTitle("nvtx",1);
143  }
144 
145 
146  //mE = dbe->book1D("E", "E", eBin, eMin, eMax);
147  //mP = dbe->book1D("P", "P", pBin, pMin, pMax);
148  // mMass = dbe->book1D("Mass", "Mass", 100, 0, 25);
149  //
150  mPhiVSEta = dbe->book2D("PhiVSEta", "PhiVSEta", 50, etaMin, etaMax, 24, phiMin, phiMax);
151  if(makedijetselection!=1){
152  mPt_1 = dbe->book1D("Pt1", "Pt1", 20, 0, 100);
153  mPt_2 = dbe->book1D("Pt2", "Pt2", 60, 0, 300);
154  mPt_3 = dbe->book1D("Pt3", "Pt3", 100, 0, 5000);
155  // Low and high pt trigger paths
156  mPt_Lo = dbe->book1D("Pt_Lo", "Pt (Pass Low Pt Jet Trigger)", 20, 0, 100);
157  //mEta_Lo = dbe->book1D("Eta_Lo", "Eta (Pass Low Pt Jet Trigger)", etaBin, etaMin, etaMax);
158  mPhi_Lo = dbe->book1D("Phi_Lo", "Phi (Pass Low Pt Jet Trigger)", phiBin, phiMin, phiMax);
159 
160  mPt_Hi = dbe->book1D("Pt_Hi", "Pt (Pass Hi Pt Jet Trigger)", 60, 0, 300);
161  mEta_Hi = dbe->book1D("Eta_Hi", "Eta (Pass Hi Pt Jet Trigger)", etaBin, etaMin, etaMax);
162  mPhi_Hi = dbe->book1D("Phi_Hi", "Phi (Pass Hi Pt Jet Trigger)", phiBin, phiMin, phiMax);
163  mNJets = dbe->book1D("NJets", "number of jets", 100, 0, 100);
164 
165  //mPt_Barrel_Lo = dbe->book1D("Pt_Barrel_Lo", "Pt Barrel (Pass Low Pt Jet Trigger)", 20, 0, 100);
166  //mPhi_Barrel_Lo = dbe->book1D("Phi_Barrel_Lo", "Phi Barrel (Pass Low Pt Jet Trigger)", phiBin, phiMin, phiMax);
167  mConstituents_Barrel = dbe->book1D("Constituents_Barrel", "Constituents Barrel above", 50, 0, 100);
168  mHFrac_Barrel = dbe->book1D("HFrac_Barrel", "HFrac Barrel", 100, 0, 1);
169  mEFrac_Barrel = dbe->book1D("EFrac_Barrel", "EFrac Barrel", 110, -0.05, 1.05);
170 
171  //mPt_EndCap_Lo = dbe->book1D("Pt_EndCap_Lo", "Pt EndCap (Pass Low Pt Jet Trigger)", 20, 0, 100);
172  //mPhi_EndCap_Lo = dbe->book1D("Phi_EndCap_Lo", "Phi EndCap (Pass Low Pt Jet Trigger)", phiBin, phiMin, phiMax);
173  mConstituents_EndCap = dbe->book1D("Constituents_EndCap", "Constituents EndCap", 50, 0, 100);
174  mHFrac_EndCap = dbe->book1D("HFrac_Endcap", "HFrac EndCap", 100, 0, 1);
175  mEFrac_EndCap = dbe->book1D("EFrac_Endcap", "EFrac EndCap", 110, -0.05, 1.05);
176 
177  //mPt_Forward_Lo = dbe->book1D("Pt_Forward_Lo", "Pt Forward (Pass Low Pt Jet Trigger)", 20, 0, 100);
178  //mPhi_Forward_Lo = dbe->book1D("Phi_Forward_Lo", "Phi Forward (Pass Low Pt Jet Trigger)", phiBin, phiMin, phiMax);
179  mConstituents_Forward = dbe->book1D("Constituents_Forward", "Constituents Forward", 50, 0, 100);
180  mHFrac_Forward = dbe->book1D("HFrac_Forward", "HFrac Forward", 100, 0, 1);
181  mEFrac_Forward = dbe->book1D("EFrac_Forward", "EFrac Forward", 110, -0.05, 1.05);
182 
183  mPt_Barrel_Hi = dbe->book1D("Pt_Barrel_Hi", "Pt Barrel (Pass Hi Pt Jet Trigger)", 60, 0, 300);
184  mPhi_Barrel_Hi = dbe->book1D("Phi_Barrel_Hi", "Phi Barrel (Pass Hi Pt Jet Trigger)", phiBin, phiMin, phiMax);
185  //mConstituents_Barrel_Hi = dbe->book1D("Constituents_Barrel_Hi", "Constituents Barrel (Pass Hi Pt Jet Trigger)", 50, 0, 100);
186  //mHFrac_Barrel_Hi = dbe->book1D("HFrac_Barrel_Hi", "HFrac Barrel (Pass Hi Pt Jet Trigger)", 100, 0, 1);
187 
188  mPt_EndCap_Hi = dbe->book1D("Pt_EndCap_Hi", "Pt EndCap (Pass Hi Pt Jet Trigger)", 60, 0, 300);
189  mPhi_EndCap_Hi = dbe->book1D("Phi_EndCap_Hi", "Phi EndCap (Pass Hi Pt Jet Trigger)", phiBin, phiMin, phiMax);
190  //mConstituents_EndCap_Hi = dbe->book1D("Constituents_EndCap_Hi", "Constituents EndCap (Pass Hi Pt Jet Trigger)", 50, 0, 100);
191  //mHFrac_EndCap_Hi = dbe->book1D("HFrac_EndCap_Hi", "HFrac EndCap (Pass Hi Pt Jet Trigger)", 100, 0, 1);
192 
193  mPt_Forward_Hi = dbe->book1D("Pt_Forward_Hi", "Pt Forward (Pass Hi Pt Jet Trigger)", 60, 0, 300);
194  mPhi_Forward_Hi = dbe->book1D("Phi_Forward_Hi", "Phi Forward (Pass Hi Pt Jet Trigger)", phiBin, phiMin, phiMax);
195  //mConstituents_Forward_Hi = dbe->book1D("Constituents_Forward_Hi", "Constituents Forward (Pass Hi Pt Jet Trigger)", 50, 0, 100);
196  //mHFrac_Forward_Hi = dbe->book1D("HFrac_Forward_Hi", "HFrac Forward (Pass Hi Pt Jet Trigger)", 100, 0, 1);
197 
198  mPhi_Barrel = dbe->book1D("Phi_Barrel", "Phi_Barrel", phiBin, phiMin, phiMax);
199  //mE_Barrel = dbe->book1D("E_Barrel", "E_Barrel", eBin, eMin, eMax);
200  mPt_Barrel = dbe->book1D("Pt_Barrel", "Pt_Barrel", ptBin, ptMin, ptMax);
201 
202  mPhi_EndCap = dbe->book1D("Phi_EndCap", "Phi_EndCap", phiBin, phiMin, phiMax);
203  //mE_EndCap = dbe->book1D("E_EndCap", "E_EndCap", eBin, eMin, 2*eMax);
204  mPt_EndCap = dbe->book1D("Pt_EndCap", "Pt_EndCap", ptBin, ptMin, ptMax);
205 
206  mPhi_Forward = dbe->book1D("Phi_Forward", "Phi_Forward", phiBin, phiMin, phiMax);
207  //mE_Forward = dbe->book1D("E_Forward", "E_Forward", eBin, eMin, 4*eMax);
208  mPt_Forward = dbe->book1D("Pt_Forward", "Pt_Forward", ptBin, ptMin, ptMax);
209 
210  // Leading Jet Parameters
211  mEtaFirst = dbe->book1D("EtaFirst", "EtaFirst", 100, -5, 5);
212  mPhiFirst = dbe->book1D("PhiFirst", "PhiFirst", 70, -3.5, 3.5);
213  //mEFirst = dbe->book1D("EFirst", "EFirst", 100, 0, 1000);
214  mPtFirst = dbe->book1D("PtFirst", "PtFirst", 100, 0, 500);
215  if(fillJIDPassFrac==1){//fillJIDPassFrac defines a collection of cleaned jets, for which we will want to fill the cleaning passing fraction
216  mLooseJIDPassFractionVSeta = dbe->bookProfile("LooseJIDPassFractionVSeta","LooseJIDPassFractionVSeta",etaBin, etaMin, etaMax,0.,1.2);
217  mLooseJIDPassFractionVSpt = dbe->bookProfile("LooseJIDPassFractionVSpt","LooseJIDPassFractionVSpt",ptBin, ptMin, ptMax,0.,1.2);
218  mTightJIDPassFractionVSeta = dbe->bookProfile("TightJIDPassFractionVSeta","TightJIDPassFractionVSeta",etaBin, etaMin, etaMax,0.,1.2);
219  mTightJIDPassFractionVSpt = dbe->bookProfile("TightJIDPassFractionVSpt","TightJIDPassFractionVSpt",ptBin, ptMin, ptMax,0.,1.2);
220 
221 
222  }
223  }
224  // CaloJet specific
225  mMaxEInEmTowers = dbe->book1D("MaxEInEmTowers", "MaxEInEmTowers", 100, 0, 100);
226  mMaxEInHadTowers = dbe->book1D("MaxEInHadTowers", "MaxEInHadTowers", 100, 0, 100);
227  if(makedijetselection!=1) {
228  mHadEnergyInHO = dbe->book1D("HadEnergyInHO", "HadEnergyInHO", 100, 0, 10);
229  mHadEnergyInHB = dbe->book1D("HadEnergyInHB", "HadEnergyInHB", 100, 0, 50);
230  mHadEnergyInHF = dbe->book1D("HadEnergyInHF", "HadEnergyInHF", 100, 0, 50);
231  mHadEnergyInHE = dbe->book1D("HadEnergyInHE", "HadEnergyInHE", 100, 0, 100);
232  mEmEnergyInEB = dbe->book1D("EmEnergyInEB", "EmEnergyInEB", 100, 0, 50);
233  mEmEnergyInEE = dbe->book1D("EmEnergyInEE", "EmEnergyInEE", 100, 0, 50);
234  mEmEnergyInHF = dbe->book1D("EmEnergyInHF", "EmEnergyInHF", 120, -20, 100);
235  }
236  mDPhi = dbe->book1D("DPhi", "dPhi btw the two leading jets", 100, 0., acos(-1.));
237 
238  //JetID variables
239 
240  mresEMF = dbe->book1D("resEMF", "resEMF", 50, 0., 1.);
241  mN90Hits = dbe->book1D("N90Hits", "N90Hits", 50, 0., 50);
242  mfHPD = dbe->book1D("fHPD", "fHPD", 50, 0., 1.);
243  mfRBX = dbe->book1D("fRBX", "fRBX", 50, 0., 1.);
244 
245  // msigmaEta = dbe->book1D("sigmaEta", "sigmaEta", 50, 0., 1.);
246  // msigmaPhi = dbe->book1D("sigmaPhi", "sigmaPhi", 50, 0., 0.5);
247 
248  if(makedijetselection==1) {
249  mDijetAsymmetry = dbe->book1D("DijetAsymmetry", "DijetAsymmetry", 100, -1., 1.);
250  mDijetBalance = dbe->book1D("DijetBalance", "DijetBalance", 100, -2., 2.);
251  if (fillJIDPassFrac==1) {
252  mLooseJIDPassFractionVSeta = dbe->bookProfile("LooseJIDPassFractionVSeta","LooseJIDPassFractionVSeta",50, -3., 3.,0.,1.2);
253  mLooseJIDPassFractionVSpt = dbe->bookProfile("LooseJIDPassFractionVSpt","LooseJIDPassFractionVSpt",ptBin, ptMin, ptMax,0.,1.2);
254  mTightJIDPassFractionVSeta = dbe->bookProfile("TightJIDPassFractionVSeta","TightJIDPassFractionVSeta",50, -3., 3.,0.,1.2);
255  mTightJIDPassFractionVSpt = dbe->bookProfile("TightJIDPassFractionVSpt","TightJIDPassFractionVSpt",ptBin, ptMin, ptMax,0.,1.2);
256  }
257  }
258 }
259 
260 
261 //void JetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup,
262 // const edm::TriggerResults& triggerResults,
263 // const reco::CaloJet& jet) {
264 
265 
266 // ***********************************************************
268  delete jetID;
269 }
270 
271 
272 // ***********************************************************
274  const reco::CaloJetCollection& caloJets,
275  const int numPV) {
276  int numofjets=0;
277  double fstPhi=0.;
278  double sndPhi=0.;
279  double diff = 0.;
280  double corr = 0.;
281  double dphi = -999. ;
282  bool thiscleaned=false;
283  bool Loosecleaned=false;
284  bool Tightcleaned=false;
285  bool thisemfclean=true;
286  bool emfcleanLoose=true;
287  bool emfcleanTight=true;
288 
289  srand( iEvent.id().event() % 10000);
290 
291 
292  if (makedijetselection == 1) {
293  //Dijet selection - careful: the pT is uncorrected!
294  //if(makedijetselection==1 && caloJets.size()>=2){
295  if(caloJets.size()>=2){
296  double dphiDJ = -999. ;
297  bool emfcleanLooseFirstJet=true;
298  bool emfcleanLooseSecondJet=true;
299  bool emfcleanTightFirstJet=true;
300  bool emfcleanTightSecondJet=true;
301  bool LoosecleanedFirstJet = false;
302  bool LoosecleanedSecondJet = false;
303  bool TightcleanedFirstJet = false;
304  bool TightcleanedSecondJet = false;
305  //both jets pass pt threshold
306  if ((caloJets.at(0)).pt() > _ptThreshold && (caloJets.at(1)).pt() > _ptThreshold ) {
307  if(fabs((caloJets.at(0)).eta())<3. && fabs((caloJets.at(1)).eta())<3. ){
308  //calculate dphi
309  dphiDJ = fabs((caloJets.at(0)).phi()-(caloJets.at(1)).phi());
310  if (dphiDJ > 3.14) dphiDJ=fabs(dphiDJ -6.28 );
311  //fill DPhi histo (before cutting)
312  if (mDPhi) mDPhi->Fill (dphiDJ);
313  //dphi cut
314  if(fabs(dphiDJ)>2.1){
315  //JetID
316  emfcleanLooseFirstJet=true;
317  emfcleanTightFirstJet=true;
318  emfcleanLooseSecondJet=true;
319  emfcleanTightSecondJet=true;
320  //jetID for first jet
321  jetID->calculate(iEvent, (caloJets.at(0)));
322  if(jetID->restrictedEMF()<_resEMFMinLoose && fabs((caloJets.at(0)).eta())<2.6) emfcleanLooseFirstJet=false;
323  if(jetID->restrictedEMF()<_resEMFMinTight && fabs((caloJets.at(0)).eta())<2.6) emfcleanTightFirstJet=false;
324  if(jetID->n90Hits()>=_n90HitsMinLoose && jetID->fHPD()<_fHPDMaxLoose && emfcleanLooseFirstJet) LoosecleanedFirstJet=true;
325  if(jetID->n90Hits()>=_n90HitsMinTight && jetID->fHPD()<_fHPDMaxTight && sqrt((caloJets.at(0)).etaetaMoment())>_sigmaEtaMinTight && sqrt((caloJets.at(0)).phiphiMoment())>_sigmaPhiMinTight && emfcleanTightFirstJet) TightcleanedFirstJet=true;
326  //fill the JID variables histograms BEFORE you cut on them
327  if (mN90Hits) mN90Hits->Fill (jetID->n90Hits());
328  if (mfHPD) mfHPD->Fill (jetID->fHPD());
329  if (mresEMF) mresEMF->Fill (jetID->restrictedEMF());
330  if (mfRBX) mfRBX->Fill (jetID->fRBX());
331 
332  //jetID for second jet
333  jetID->calculate(iEvent, (caloJets.at(1)));
334  if(jetID->restrictedEMF()<_resEMFMinLoose && fabs((caloJets.at(1)).eta())<2.6) emfcleanLooseSecondJet=false;
335  if(jetID->restrictedEMF()<_resEMFMinTight && fabs((caloJets.at(1)).eta())<2.6) emfcleanTightSecondJet=false;
336  if(jetID->n90Hits()>=_n90HitsMinLoose && jetID->fHPD()<_fHPDMaxLoose && emfcleanLooseSecondJet) LoosecleanedSecondJet=true;
337  if(jetID->n90Hits()>=_n90HitsMinTight && jetID->fHPD()<_fHPDMaxTight && sqrt((caloJets.at(1)).etaetaMoment())>_sigmaEtaMinTight && sqrt((caloJets.at(1)).phiphiMoment())>_sigmaPhiMinTight && emfcleanTightSecondJet) TightcleanedSecondJet=true;
338  //fill the JID variables histograms BEFORE you cut on them
339  if (mN90Hits) mN90Hits->Fill (jetID->n90Hits());
340  if (mfHPD) mfHPD->Fill (jetID->fHPD());
341  if (mresEMF) mresEMF->Fill (jetID->restrictedEMF());
342  if (mfRBX) mfRBX->Fill (jetID->fRBX());
343 
344  if(LoosecleanedFirstJet && LoosecleanedSecondJet) { //only if both jets are (loose) cleaned
345  //fill histos for first jet
346  if (mPt) mPt->Fill ((caloJets.at(0)).pt());
347  if (mEta) mEta->Fill ((caloJets.at(0)).eta());
348  if (mPhi) mPhi->Fill ((caloJets.at(0)).phi());
349  if (mPhiVSEta) mPhiVSEta->Fill((caloJets.at(0)).eta(),(caloJets.at(0)).phi());
350  if (mConstituents) mConstituents->Fill ((caloJets.at(0)).nConstituents());
351  if (mHFrac) mHFrac->Fill ((caloJets.at(0)).energyFractionHadronic());
352  if (mEFrac) mEFrac->Fill ((caloJets.at(0)).emEnergyFraction());
353  //if (mE) mE->Fill ((caloJets.at(0)).energy());
354  //if (mP) mP->Fill ((caloJets.at(0)).p());
355  //if (mMass) mMass->Fill ((caloJets.at(0)).mass());
356  if (mMaxEInEmTowers) mMaxEInEmTowers->Fill ((caloJets.at(0)).maxEInEmTowers());
357  if (mMaxEInHadTowers) mMaxEInHadTowers->Fill ((caloJets.at(0)).maxEInHadTowers());
358  if (mN90Hits) mN90Hits->Fill (jetID->n90Hits());
359  if (mfHPD) mfHPD->Fill (jetID->fHPD());
360  if (mresEMF) mresEMF->Fill (jetID->restrictedEMF());
361  if (mfRBX) mfRBX->Fill (jetID->fRBX());
362  //sigmaeta and sigmaphi only used in the tight selection.
363  //fill the histos for them AFTER the loose selection
364  // if (msigmaEta) msigmaEta->Fill(sqrt((caloJets.at(0)).etaetaMoment()));
365  // if (msigmaPhi) msigmaPhi->Fill(sqrt((caloJets.at(0)).phiphiMoment()));
366  //fill histos for second jet
367  if (mPt) mPt->Fill ((caloJets.at(1)).pt());
368  if (mEta) mEta->Fill ((caloJets.at(1)).eta());
369  if (mPhi) mPhi->Fill ((caloJets.at(1)).phi());
370  if (mPhiVSEta) mPhiVSEta->Fill((caloJets.at(1)).eta(),(caloJets.at(1)).phi());
371  if (mConstituents) mConstituents->Fill ((caloJets.at(1)).nConstituents());
372  if (mHFrac) mHFrac->Fill ((caloJets.at(1)).energyFractionHadronic());
373  if (mEFrac) mEFrac->Fill ((caloJets.at(1)).emEnergyFraction());
374  //if (mE) mE->Fill ((caloJets.at(1)).energy());
375  //if (mP) mP->Fill ((caloJets.at(1)).p());
376  //if (mMass) mMass->Fill ((caloJets.at(1)).mass());
377  if (mMaxEInEmTowers) mMaxEInEmTowers->Fill ((caloJets.at(1)).maxEInEmTowers());
378  if (mMaxEInHadTowers) mMaxEInHadTowers->Fill ((caloJets.at(1)).maxEInHadTowers());
379  //sigmaeta and sigmaphi only used in the tight selection.
380  //fill the histos for them AFTER the loose selection
381  // if (msigmaEta) msigmaEta->Fill(sqrt((caloJets.at(1)).etaetaMoment()));
382  // if (msigmaPhi) msigmaPhi->Fill(sqrt((caloJets.at(1)).phiphiMoment()));
383 
384 
385  // Fill NPV profiles
386  //----------------------------------------------------------------
387  for (int ijet=0; ijet<2; ijet++) {
388 
389  if (mPt_profile) mPt_profile ->Fill(numPV, (caloJets.at(ijet)).pt());
390  if (mEta_profile) mEta_profile ->Fill(numPV, (caloJets.at(ijet)).eta());
391  if (mPhi_profile) mPhi_profile ->Fill(numPV, (caloJets.at(ijet)).phi());
392  if (mConstituents_profile) mConstituents_profile->Fill(numPV, (caloJets.at(ijet)).nConstituents());
393  if (mHFrac_profile) mHFrac_profile ->Fill(numPV, (caloJets.at(ijet)).energyFractionHadronic());
394  if (mEFrac_profile) mEFrac_profile ->Fill(numPV, (caloJets.at(ijet)).emEnergyFraction());
395  }
396  }
397 
398 
399  //let's see how many of these jets passed the JetID cleaning
400  if(fillJIDPassFrac==1) {
401  if(LoosecleanedFirstJet) {
402  mLooseJIDPassFractionVSeta->Fill((caloJets.at(0)).eta(),1.);
403  mLooseJIDPassFractionVSpt->Fill((caloJets.at(0)).pt(),1.);
404  } else {
405  mLooseJIDPassFractionVSeta->Fill((caloJets.at(0)).eta(),0.);
406  mLooseJIDPassFractionVSpt->Fill((caloJets.at(0)).pt(),0.);
407  }
408  if(LoosecleanedSecondJet) {
409  mLooseJIDPassFractionVSeta->Fill((caloJets.at(1)).eta(),1.);
410  mLooseJIDPassFractionVSpt->Fill((caloJets.at(1)).pt(),1.);
411  } else {
412  mLooseJIDPassFractionVSeta->Fill((caloJets.at(1)).eta(),0.);
413  mLooseJIDPassFractionVSpt->Fill((caloJets.at(1)).pt(),0.);
414  }
415  //TIGHT JID
416  if(TightcleanedFirstJet) {
417  mTightJIDPassFractionVSeta->Fill((caloJets.at(0)).eta(),1.);
418  mTightJIDPassFractionVSpt->Fill((caloJets.at(0)).pt(),1.);
419  } else {
420  mTightJIDPassFractionVSeta->Fill((caloJets.at(0)).eta(),0.);
421  mTightJIDPassFractionVSpt->Fill((caloJets.at(0)).pt(),0.);
422  }
423  if(TightcleanedSecondJet) {
424  mTightJIDPassFractionVSeta->Fill((caloJets.at(1)).eta(),1.);
425  mTightJIDPassFractionVSpt->Fill((caloJets.at(1)).pt(),1.);
426  } else {
427  mTightJIDPassFractionVSeta->Fill((caloJets.at(1)).eta(),0.);
428  mTightJIDPassFractionVSpt->Fill((caloJets.at(1)).pt(),0.);
429  }
430 
431  }//if fillJIDPassFrac
432  }// FABS DPHI < 2.1
433  }// fabs eta < 3
434  }// pt jets > threshold
435  //now do the dijet balance and asymmetry calculations
436  if (fabs(caloJets.at(0).eta() < 1.4)) {
437  double pt_dijet = (caloJets.at(0).pt() + caloJets.at(1).pt())/2;
438 
439  double dPhi = fabs((caloJets.at(0)).phi()-(caloJets.at(1)).phi());
440  if (dPhi > 3.14) dPhi=fabs(dPhi -6.28 );
441 
442  if (dPhi > 2.7) {
443  double pt_probe;
444  double pt_barrel;
445  int jet1, jet2;
446 
447  int randJet = rand() % 2;
448 
449  if (fabs(caloJets.at(1).eta() < 1.4)) {
450  if (randJet) {
451  jet1 = 0;
452  jet2 = 1;
453  }
454  else {
455  jet1 = 1;
456  jet2 = 0;
457  }
458 
459  /***Di-Jet Asymmetry****
460  * leading jets eta < 1.4
461  * leading jets dphi > 2.7
462  * pt_third jet < threshold
463  * A = (pt_1 - pt_2)/(pt_1 + pt_2)
464  * jets 1 and two are randomly ordered
465  */
466  bool thirdJetCut = true;
467  for (unsigned int third = 2; third < caloJets.size(); ++third)
468  if (caloJets.at(third).pt() > _asymmetryThirdJetCut)
469  thirdJetCut = false;
470  if (thirdJetCut) {
471  double dijetAsymmetry = (caloJets.at(jet1).pt() - caloJets.at(jet2).pt()) / (caloJets.at(jet1).pt() + caloJets.at(jet2).pt());
472  mDijetAsymmetry->Fill(dijetAsymmetry);
473  }// end restriction on third jet pt in asymmetry calculation
474 
475  }
476  else {
477  jet1 = 0;
478  jet2 = 1;
479  }
480 
481  pt_barrel = caloJets.at(jet1).pt();
482  pt_probe = caloJets.at(jet2).pt();
483 
484  //dijet balance cuts
485  /***Di-Jet Balance****
486  * pt_dijet = (pt_probe+pt_barrel)/2
487  * leading jets dphi > 2.7
488  * reject evnets where pt_third/pt_dijet > 0.2
489  * pv selection
490  * B = (pt_probe - pt_barrel)/pt_dijet
491  * select probe randomly from 2 jets if both leading jets are in the barrel
492  */
493  bool thirdJetCut = true;
494  for (unsigned int third = 2; third < caloJets.size(); ++third)
495  if (caloJets.at(third).pt()/pt_dijet > _balanceThirdJetCut)
496  thirdJetCut = false;
497  if (thirdJetCut) {
498  double dijetBalance = (pt_probe - pt_barrel) / pt_dijet;
499  mDijetBalance->Fill(dijetBalance);
500  }// end restriction on third jet pt ratio in balance calculation
501  }// dPhi > 2.7
502  }// leading jet eta cut for asymmetry and balance calculations
503  }//jet size >= 2
504  }// do dijet selection
505  else {
506  for (reco::CaloJetCollection::const_iterator jet = caloJets.begin(); jet!=caloJets.end(); ++jet) {
507  LogTrace(jetname)<<"[JetAnalyzer] Analyze Calo Jet";
508  Loosecleaned=false;
509  Tightcleaned=false;
510  if (jet == caloJets.begin()) {
511  fstPhi = jet->phi();
512  _leadJetFlag = 1;
513  } else {
514  _leadJetFlag = 0;
515  }
516  if (jet == (caloJets.begin()+1)) sndPhi = jet->phi();
517  //jetID
518  jetID->calculate(iEvent, *jet);
519  //minimal (uncorrected!) pT cut
520  if (jet->pt() > _ptThreshold) {
521  // if (msigmaEta) msigmaEta->Fill(sqrt(jet->etaetaMoment()));
522  // if (msigmaPhi) msigmaPhi->Fill(sqrt(jet->phiphiMoment()));
523  //cleaning to use for filling histograms
524  thisemfclean=true;
525  if(jetID->restrictedEMF()<_resEMFMin && fabs(jet->eta())<2.6) thisemfclean=false;
526  if(jetID->n90Hits()>=_n90HitsMin && jetID->fHPD()<_fHPDMax && thisemfclean) thiscleaned=true;
527  //loose and tight cleaning, used to fill the JetIDPAssFraction histos
528  if(jetID->n90Hits()>=_n90HitsMinLoose && jetID->fHPD()<_fHPDMaxLoose && emfcleanLoose) Loosecleaned=true;
529  if(jetID->n90Hits()>=_n90HitsMinTight && jetID->fHPD()<_fHPDMaxTight && sqrt(jet->etaetaMoment())>_sigmaEtaMinTight && sqrt(jet->phiphiMoment())>_sigmaPhiMinTight && emfcleanTight) Tightcleaned=true;
530 
531  if(fillJIDPassFrac==1) {
532  if(Loosecleaned) {
533  mLooseJIDPassFractionVSeta->Fill(jet->eta(),1.);
534  mLooseJIDPassFractionVSpt->Fill(jet->pt(),1.);
535  } else {
536  mLooseJIDPassFractionVSeta->Fill(jet->eta(),0.);
537  mLooseJIDPassFractionVSpt->Fill(jet->pt(),0.);
538  }
539  //TIGHT
540  if(Tightcleaned) {
541  mTightJIDPassFractionVSeta->Fill(jet->eta(),1.);
542  mTightJIDPassFractionVSpt->Fill(jet->pt(),1.);
543  } else {
544  mTightJIDPassFractionVSeta->Fill(jet->eta(),0.);
545  mTightJIDPassFractionVSpt->Fill(jet->pt(),0.);
546  }
547  }
548  //eventually we could define the "cleaned" flag differently for e.g. HF
549  if(thiscleaned) {
550  numofjets++ ;
551  jetME->Fill(1);
552 
553  // Leading jet
554  // Histograms are filled once per event
555  if (_leadJetFlag == 1) {
556  if (mEtaFirst) mEtaFirst->Fill (jet->eta());
557  if (mPhiFirst) mPhiFirst->Fill (jet->phi());
558  //if (mEFirst) mEFirst->Fill (jet->energy());
559  if (mPtFirst) mPtFirst->Fill (jet->pt());
560  }
561  // --- Passed the low pt jet trigger
562  if (_JetLoPass == 1) {
563  /* if (fabs(jet->eta()) <= 1.3) {
564  if (mPt_Barrel_Lo) mPt_Barrel_Lo->Fill(jet->pt());
565  if (mEta_Lo) mEta_Lo->Fill(jet->eta());
566  if (mPhi_Barrel_Lo) mPhi_Barrel_Lo->Fill(jet->phi());
567  }
568  if ( (fabs(jet->eta()) > 1.3) && (fabs(jet->eta()) <= 3) ) {
569  if (mPt_EndCap_Lo) mPt_EndCap_Lo->Fill(jet->pt());
570  if (mEta_Lo) mEta_Lo->Fill(jet->eta());
571  if (mPhi_EndCap_Lo) mPhi_EndCap_Lo->Fill(jet->phi());
572  }
573  if (fabs(jet->eta()) > 3.0) {
574  if (mPt_Forward_Lo) mPt_Forward_Lo->Fill(jet->pt());
575  if (mEta_Lo) mEta_Lo->Fill(jet->eta());
576  if (mPhi_Forward_Lo) mPhi_Forward_Lo->Fill(jet->phi());
577  } */
578  //if (mEta_Lo) mEta_Lo->Fill (jet->eta());
579  if (mPhi_Lo) mPhi_Lo->Fill (jet->phi());
580  if (mPt_Lo) mPt_Lo->Fill (jet->pt());
581  }
582 
583  // --- Passed the high pt jet trigger
584  if (_JetHiPass == 1) {
585  if (fabs(jet->eta()) <= 1.3) {
586  if (mPt_Barrel_Hi && jet->pt()>100.) mPt_Barrel_Hi->Fill(jet->pt());
587  if (mEta_Hi && jet->pt()>100.) mEta_Hi->Fill(jet->eta());
588  if (mPhi_Barrel_Hi) mPhi_Barrel_Hi->Fill(jet->phi());
589  //if (mConstituents_Barrel_Hi) mConstituents_Barrel_Hi->Fill(jet->nConstituents());
590  //if (mHFrac_Barrel_Hi) mHFrac_Barrel_Hi->Fill(jet->energyFractionHadronic());
591  }
592  if ( (fabs(jet->eta()) > 1.3) && (fabs(jet->eta()) <= 3) ) {
593  if (mPt_EndCap_Hi && jet->pt()>100.) mPt_EndCap_Hi->Fill(jet->pt());
594  if (mEta_Hi && jet->pt()>100.) mEta_Hi->Fill(jet->eta());
595  if (mPhi_EndCap_Hi) mPhi_EndCap_Hi->Fill(jet->phi());
596  //if (mConstituents_EndCap_Hi) mConstituents_EndCap_Hi->Fill(jet->nConstituents());
597  //if (mHFrac_EndCap_Hi) mHFrac_EndCap_Hi->Fill(jet->energyFractionHadronic());
598  }
599  if (fabs(jet->eta()) > 3.0) {
600  if (mPt_Forward_Hi && jet->pt()>100.) mPt_Forward_Hi->Fill(jet->pt());
601  if (mEta_Hi && jet->pt()>100.) mEta_Hi->Fill(jet->eta());
602  if (mPhi_Forward_Hi) mPhi_Forward_Hi->Fill(jet->phi());
603  //if (mConstituents_Forward_Hi) mConstituents_Forward_Hi->Fill(jet->nConstituents());
604  //if (mHFrac_Forward_Hi) mHFrac_Forward_Hi->Fill(jet->energyFractionHadronic());
605  }
606 
607  if (mEta_Hi && jet->pt()>100.) mEta_Hi->Fill (jet->eta());
608  if (mPhi_Hi) mPhi_Hi->Fill (jet->phi());
609  if (mPt_Hi) mPt_Hi->Fill (jet->pt());
610  }
611 
612  if (mPt) mPt->Fill (jet->pt());
613  if (mPt_1) mPt_1->Fill (jet->pt());
614  if (mPt_2) mPt_2->Fill (jet->pt());
615  if (mPt_3) mPt_3->Fill (jet->pt());
616  if (mEta) mEta->Fill (jet->eta());
617  if (mPhi) mPhi->Fill (jet->phi());
618 
619  if (mPhiVSEta) mPhiVSEta->Fill(jet->eta(),jet->phi());
620 
621  if (mConstituents) mConstituents->Fill (jet->nConstituents());
622  if (mHFrac) mHFrac->Fill (jet->energyFractionHadronic());
623  if (mEFrac) mEFrac->Fill (jet->emEnergyFraction());
624 
625 
626  // Fill NPV profiles
627  //--------------------------------------------------------------------
628  if (mPt_profile) mPt_profile ->Fill(numPV, jet->pt());
629  if (mEta_profile) mEta_profile ->Fill(numPV, jet->eta());
630  if (mPhi_profile) mPhi_profile ->Fill(numPV, jet->phi());
631  if (mConstituents_profile) mConstituents_profile->Fill(numPV, jet->nConstituents());
632  if (mHFrac_profile) mHFrac_profile ->Fill(numPV, jet->energyFractionHadronic());
633  if (mEFrac_profile) mEFrac_profile ->Fill(numPV, jet->emEnergyFraction());
634 
635 
636  if (fabs(jet->eta()) <= 1.3) {
637  if (mPt_Barrel) mPt_Barrel->Fill (jet->pt());
638  if (mPhi_Barrel) mPhi_Barrel->Fill (jet->phi());
639  //if (mE_Barrel) mE_Barrel->Fill (jet->energy());
640  if (mConstituents_Barrel) mConstituents_Barrel->Fill(jet->nConstituents());
641  if (mHFrac_Barrel) mHFrac_Barrel->Fill(jet->energyFractionHadronic());
642  if (mEFrac_Barrel) mEFrac_Barrel->Fill(jet->emEnergyFraction());
643  }
644  if ( (fabs(jet->eta()) > 1.3) && (fabs(jet->eta()) <= 3) ) {
645  if (mPt_EndCap) mPt_EndCap->Fill (jet->pt());
646  if (mPhi_EndCap) mPhi_EndCap->Fill (jet->phi());
647  //if (mE_EndCap) mE_EndCap->Fill (jet->energy());
648  if (mConstituents_EndCap) mConstituents_EndCap->Fill(jet->nConstituents());
649  if (mHFrac_EndCap) mHFrac_EndCap->Fill(jet->energyFractionHadronic());
650  if (mEFrac_EndCap) mEFrac_EndCap->Fill(jet->emEnergyFraction());
651  }
652  if (fabs(jet->eta()) > 3.0) {
653  if (mPt_Forward) mPt_Forward->Fill (jet->pt());
654  if (mPhi_Forward) mPhi_Forward->Fill (jet->phi());
655  //if (mE_Forward) mE_Forward->Fill (jet->energy());
656  if (mConstituents_Forward) mConstituents_Forward->Fill(jet->nConstituents());
657  if (mHFrac_Forward) mHFrac_Forward->Fill(jet->energyFractionHadronic());
658  if (mEFrac_Forward) mEFrac_Forward->Fill(jet->emEnergyFraction());
659  }
660 
661  //if (mE) mE->Fill (jet->energy());
662  //if (mP) mP->Fill (jet->p());
663  // if (mMass) mMass->Fill (jet->mass());
664 
665  if (mMaxEInEmTowers) mMaxEInEmTowers->Fill (jet->maxEInEmTowers());
666  if (mMaxEInHadTowers) mMaxEInHadTowers->Fill (jet->maxEInHadTowers());
667 
668  if (mHadEnergyInHO) mHadEnergyInHO->Fill (jet->hadEnergyInHO());
669  if (mHadEnergyInHB) mHadEnergyInHB->Fill (jet->hadEnergyInHB());
670  if (mHadEnergyInHF) mHadEnergyInHF->Fill (jet->hadEnergyInHF());
671  if (mHadEnergyInHE) mHadEnergyInHE->Fill (jet->hadEnergyInHE());
672  if (mEmEnergyInEB) mEmEnergyInEB->Fill (jet->emEnergyInEB());
673  if (mEmEnergyInEE) mEmEnergyInEE->Fill (jet->emEnergyInEE());
674  if (mEmEnergyInHF) mEmEnergyInHF->Fill (jet->emEnergyInHF());
675 
676  if (mN90Hits) mN90Hits->Fill (jetID->n90Hits());
677  if (mfHPD) mfHPD->Fill (jetID->fHPD());
678  if (mresEMF) mresEMF->Fill (jetID->restrictedEMF());
679  if (mfRBX) mfRBX->Fill (jetID->fRBX());
680 
681  //calculate correctly the dphi
682  if(numofjets>1) {
683  diff = fabs(fstPhi - sndPhi);
684  corr = 2*acos(-1.) - diff;
685  if(diff < acos(-1.)) {
686  dphi = diff;
687  } else {
688  dphi = corr;
689  }
690  }
691  }
692  }//pt cut
693  }
694 
695  if (mNJets) mNJets->Fill (numofjets);
696  if (mDPhi && dphi>-998.) mDPhi->Fill (dphi);
697 
698 
699  if (mNJets_profile) mNJets_profile->Fill(numPV, numofjets);
700 
701 
702  }//not dijet
703 }
704 
EventNumber_t event() const
Definition: EventID.h:44
dictionary parameters
Definition: Parameters.py:2
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:717
JetAnalyzer(const edm::ParameterSet &)
Constructor.
Definition: JetAnalyzer.cc:22
void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
T eta() const
const int nbinsPV
const double PVup
int iEvent
Definition: GenABIO.cc:243
virtual ~JetAnalyzer()
Destructor.
Definition: JetAnalyzer.cc:46
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
T sqrt(T t)
Definition: SSEVec.h:46
MonitorElement * bookProfile(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, const char *option="s")
Definition: DQMStore.cc:1031
#define LogTrace(id)
const double PVlow
JetCorrectorParameters corr
Definition: classes.h:9
void endJob()
Finish up a job.
Definition: JetAnalyzer.cc:267
void beginJob(DQMStore *dbe)
Inizialize parameters for histo binning.
Definition: JetAnalyzer.cc:50
edm::EventID id() const
Definition: EventBase.h:56
Signal rand(Signal arg)
Definition: vlib.cc:442
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:845
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:429
void analyze(const edm::Event &, const edm::EventSetup &, const reco::CaloJetCollection &caloJets, const int numPV)
Get the analysis.
Definition: JetAnalyzer.cc:273
std::vector< CaloJet > CaloJetCollection
collection of CaloJet objects
Definition: DDAxes.h:10