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