ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SCerci/src/ForwardJet_coneJ.cc
Revision: 1.2
Committed: Fri Oct 5 08:01:48 2007 UTC (17 years, 7 months ago) by salim
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -1 lines
Error occurred while calculating annotation data.
Log Message:
*** empty log message ***

File Contents

# Content
1 #include <iostream>
2 #include "Analysis/ForwardJet/interface/ForwardJet.h"
3 #include "DataFormats/Common/interface/Ref.h"
4 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
5 #include "DataFormats/Common/interface/Ref.h"
6 #include "DataFormats/JetReco/interface/Jet.h"
7 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
8 #include "DataFormats/JetReco/interface/BasicJetCollection.h"
9 #include "DataFormats/JetReco/interface/GenJet.h"
10 #include "DataFormats/JetReco/interface/CaloJet.h"
11 #include "DataFormats/METReco/interface/CaloMET.h"
12 #include "DataFormats/METReco/interface/GenMET.h"
13 #include "DataFormats/JetReco/interface/GenJetfwd.h"
14 #include "DataFormats/JetReco/interface/BasicJet.h"
15 #include "DataFormats/JetReco/interface/BasicJetfwd.h"
16 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
17 #include "CLHEP/HepMC/ReadHepMC.h"
18 #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
19 #include "FWCore/Framework/interface/Handle.h"
20 #include "FWCore/Framework/interface/Event.h"
21 #include "FWCore/ParameterSet/interface/ParameterSet.h"
22 //Jet_correct
23 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
24 //Jet_id
25 #include "RecoBTag/MCTools/interface/MCParticleInfo.h"
26 #include "RecoBTag/MCTools/interface/MCBaseParticle.h"
27 #include "RecoBTag/MCTools/interface/MCParton.h"
28 #include "DataFormats/BTauReco/interface/JetTag.h"
29 #include "RecoBTag/MCTools/interface/JetFlavourIdentifier.h"
30 //
31
32
33 #include "TFile.h"
34 #include "TTree.h"
35 #include "TH1.h"
36 #include "TH1F.h"
37 #include "TH2F.h"
38 #include "TMath.h"
39 #include "TF1.h"
40 #include "TCanvas.h"
41 #include "TGraphErrors.h"
42 #include "TProfile.h"
43
44 using namespace edm;
45 using namespace std;
46 using namespace reco;
47
48
49 int eventnum=0;
50 class GenJetSort{
51 public:
52 bool operator()(const GenJet& a, const GenJet& b) {
53 return a.et() > b.et();
54 }
55 };
56
57 class CaloJetSort{
58 public:
59 bool operator()(const CaloJet& a, const CaloJet& b) {
60 return a.et() > b.et();
61 }
62 };
63 //class CorJetSort{
64 //public:
65 //bool operator()(const CorJet& a, const CorJet& b) {
66 // return a.et() > b.et();
67 //}
68 //};
69 //ForwardJet::ForwardJet( const ParameterSet& pset )
70 //{}
71
72 ForwardJet::ForwardJet( const ParameterSet& pset ) :
73 //fOutputFileName( pset.getUntrackedParameter<string>("HistOutFile",string("OutPutFile"))),
74 fOutputFileName( pset.getParameter<string>("OutPutFile")),
75 //fOutputFileName( pset.getParameter<string>("HistOutFile",string("/data2/salim/ForwardJet/ForwardJet_QCD_pt_120_170.root"))),
76 fOutputFile(0),
77
78 CorJetAlgorithm( pset.getParameter<string>( "CorJetAlgorithm" ) )
79
80
81 {
82 jfi = JetFlavourIdentifier(pset.getParameter<edm::ParameterSet>("jetIdParameters"));
83 }
84
85 void ForwardJet::beginJob( const EventSetup& )
86 {
87 fOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
88
89 //RecGenTree = new TTree("RecJet_GenJet", "Analysis");
90 //RecGenTree->Branch("fNRecJ", &fNRecJ, "fNRecJ/I");
91 //RecGenTree->Branch("fRecJetPt", &fRecJetPt, "fRecJetPt[fNRecJ]/D");
92
93 // limits for histograms
94 int nbinspt = 250;
95 int nbinsphi = 100;
96 int nbinsE= 250;
97 int nbinsEta =100;
98 int nbinsBjoX=600;
99 double ptmin = 0; //GeV
100 double ptmax = 250; //GeV
101 double ymin = -10;
102 double ymax = 10;
103 double BjoXmin = -6;
104 double BjoXmax = 0;
105 double phimin = -TMath::Pi();
106 double phimax = TMath::Pi();
107 double Emin=0; //GeV
108 double Emax=250; //GeV
109 double ETCut=20; //GeV
110 double threshold = 0.05;
111 double threshold_1 = 0.04;// min distance between parton and GenJEt
112 double sqrts = 14000.;//Gev
113 float BjoX;//Bjorken-x
114 float Min_deltaR;// minumum distance between parton and GenJEt
115 // create histograms
116 ptGen = new TH1D("ptGen","QCD GenJet pT Pt_15_120 Gev",nbinspt,ptmin,ptmax);
117 ptGen->Sumw2();
118 ptParton = new TH1D("ptPArton","QCD Parton pT ",nbinspt,ptmin,ptmax);
119 ptParton->Sumw2();
120 yGen = new TH1D("yGen","QCD GenJet rapidity Pt_15_120 Gev ",nbinsEta,ymin,ymax);
121 yGen->Sumw2();
122 phiGen = new TH1D("phiGen","QCD GenJet Phi Pt_15_120 Gev",nbinsphi,phimin,phimax);
123 phiGen->Sumw2();
124 EGen = new TH1D ("EGen","QCD GenJet E Pt_15_120 Gev ",nbinsE,Emin, Emax);
125 EGen->Sumw2();
126
127 ptCalo = new TH1D("ptCalo","QCD Calo Pt Pt_15_120 Gev",nbinspt,ptmin,ptmax);
128 ptCalo->Sumw2();
129
130 corr_ptCalo = new TH1D("corr_ptCalo","Corrected Calo QCD jets pT Pt_15_120 Gev ",nbinspt,ptmin,ptmax);
131 corr_ptCalo->Sumw2();
132
133 yCalo = new TH1D("yCalo","Calo QCD jets rapidity Pt_15_120 Gev",nbinsEta,ymin,ymax);
134 yCalo->Sumw2();
135 phiCalo = new TH1D("phiCalo","Calo QCD jets phi Pt_15_120 Gev ",nbinsphi,phimin,phimax);
136 phiCalo->Sumw2();
137 ECalo = new TH1D ("ECalo","Calo QCD jets energy Pt_15_120 Gev ",nbinsE,Emin, Emax);
138 ECalo->Sumw2();
139 corr_ECalo = new TH1D ("corr_ECalo","Corrected Calo QCD jets Pt_15_120 Gev ",nbinsE,Emin, Emax);
140 corr_ECalo->Sumw2();
141
142 trueRecoJet = new TH1D ("trueRecoJet","True Reco Jet for E_{T}>20 GeV ",nbinsE,Emin, Emax);
143 trueRecoJet->Sumw2();
144
145 RatioE_vs_Eta= new TProfile(" RatioE_vs_Eta","E^{reco}_{T}/E^{MC}_{T} for QCD jets Pt_15_120 Gev ", (nbinsEta/2), 0.,5.);
146
147 RatioE_vs_E= new TProfile("RatioE_vs_E ","E^{reco}_{T}/E^{MC}_{T} for QCD jets Pt_15_120 Gev", nbinsE,Emin ,Emax);
148
149 //HF
150 pid_highx_parton_01_hf = new TH1D("pid_highx_parton01_hf","PID of parton (jet)with x>0.1 for HF",30,-0.5,29.5);
151 pid_highx_parton_02_hf = new TH1D("pid_highx_parton02_hf","PID of parton (jet)with x>0.2 for HF",30,-0.5,29.5);
152 pid_highx_parton_04_hf = new TH1D("pid_highx_parton04_hf","PID of parton (jet)with x>0.4 for HF",30,-0.5,29.5);
153 pid_highx_parton_06_hf = new TH1D("pid_highx_parton06_hf","PID of parton (jet)with x>0.6 for HF",30,-0.5,29.5);
154 pid_Allx_parton_hf= new TH1D("pid_Allx_parton_hf","PID of parton (jet) 10^{-5}<x<0 for HF",30,-0.5,29.5);
155 //CASTOR
156 pid_highx_parton_01_castor = new TH1D("pid_highx_parton01_castor","PID of parton (jet)with x>0.1 for CASTOR",30,-0.5,29.5);
157 pid_highx_parton_02_castor = new TH1D("pid_highx_parton02_castor","PID of parton (jet)with x>0.2 for CASTOR",30,-0.5,29.5);
158 pid_highx_parton_04_castor = new TH1D("pid_highx_parton04_castor","PID of parton (jet)with x>0.4 for CASTOR",30,-0.5,29.5);
159 pid_highx_parton_06_castor = new TH1D("pid_highx_parton06_castor","PID of parton (jet)with x>0.6 for CASTOR",30,-0.5,29.5);
160 pid_Allx_parton_CASTOR= new TH1D("pid_Allx_parton_CASTOR","PID of parton (jet) 10^{-6}<x<0 for CASTOR",30,-0.5,29.5);
161 //DeltaR & GenJetEt-PartonEt
162 Min_DeltaR=new TH1D("Min_DeltaR","Min_deltaR for Parton and GenJet ",100,-0.01,1.);
163 Min_DeltaR_02=new TH2D("Min_DeltaR_02","GenJetP_{T}-PartonP_{T} for #DeltaR=0.02",nbinsE,Emin,Emax,400,-200.,200.);
164 Min_DeltaR_04=new TH2D("Min_DeltaR_04","GenJetP_{T}-PartonP_{T} for #DeltaR=0.04",nbinsE,Emin,Emax,400,-200.,200.);
165 Min_DeltaR_06=new TH2D("Min_DeltaR_06","GenJetP_{T}-PartonP_{T} for #DeltaR=0.06",nbinsE,Emin,Emax,400,-200.,200.);
166 Min_DeltaR_08=new TH2D("Min_DeltaR_08","GenJetP_{T}-PartonP_{T} for #DeltaR=0.08",nbinsE,Emin,Emax,400,-200.,200.);
167 Min_DeltaR_1=new TH2D("Min_DeltaR_1","GenJetP_{T}-PartonP_{T} for #DeltaR=0.1",nbinsE,Emin,Emax,400,-200.,200.);
168 Min_DeltaR_2=new TH2D("Min_DeltaR_2","GenJetP_{T}-PartonP_{T} for #DeltaR=0.2",nbinsE,Emin,Emax,400,-200.,200.);
169 Min_DeltaR_5=new TH2D("Min_DeltaR_5","GenJetP_{T}-PartonP_{T} for #DeltaR=0.5",nbinsE,Emin,Emax,400,-200.,200.);
170 Min_DeltaR_8=new TH2D("Min_DeltaR_8","GenJetP_{T}-PartonP_{T} for #DeltaR=0.8",nbinsE,Emin,Emax,400,-200.,200.);
171 //Creating a binning in ET and an histogram for each of the ET bin
172 unsigned int nbreBins=11;
173 //EThist [11];
174 char nameHist [50];//for Et distribution
175 char nameHist1 [50];//for reco jet phi distribution
176 char nameHist2 [50];//for reco jet eta distribution
177 char label [500];//for Et distribution
178 char label1 [500];//for reco jet phi distribution
179 char label2 [500];//for reco jet eta distribution
180
181 double limBins[] = {10,20,30,40,50,60,70,80,100.0,120.,150.,250};
182
183 ratioE = new TH1D ("ratioE", "E^{reco}_{T}/E^{MC}_{T}",nbreBins, Emin, Emax);
184 ratioE->Sumw2();
185 ratioE->SetBins(nbreBins,limBins);
186
187 sigmaErEm = new TH1D ("sigmaErEm", "sigma Ereco / Emc",nbreBins, Emin, Emax);
188
189 sigmaErEm->Sumw2();
190 sigmaErEm->SetBins(nbreBins,limBins);
191
192 for (unsigned int indexH=0;indexH<nbreBins;indexH++)
193 {
194 sprintf(nameHist,"ETdistribution%i",indexH);
195 sprintf(nameHist1,"phidistribution%i",indexH);
196 sprintf(nameHist2,"etadistribution%i",indexH);
197
198 sprintf(label,"E_{T}(reco)/E_{T}(Gen) #in #left[%f , %f #right] GeV",limBins[indexH],limBins[indexH+1]);
199 sprintf(label1,"#phi_{rec} #in #left[%f , %f #right] GeV",limBins[indexH],limBins[indexH+1]);
200 sprintf(label2,"#eta_{rec} #in #left[%f , %f #right] GeV",limBins[indexH],limBins[indexH+1]);
201
202 EThist[indexH]=new TH1D(nameHist,label,nbinsE,0.2, 1.8);
203 EThist[indexH]->Sumw2();
204 EThist[indexH]->GetXaxis()->SetTitle("E_{T}^{reco}/E_{t}^{gen}");
205 EThist[indexH]->GetYaxis()->SetTitle("#");
206
207 Phihist[indexH]=new TH1D(nameHist1,label1,480, 0.2, 1.8 );
208 Phihist[indexH]->Sumw2();
209 Phihist[indexH]->GetXaxis()->SetTitle("#phi_{rec}/#phi_{gen}");
210 Phihist[indexH]->GetYaxis()->SetTitle("#");
211
212 Etahist[indexH]=new TH1D(nameHist2,label2,400,0.8, 1.2);
213 Etahist[indexH]->Sumw2();
214 Etahist[indexH]->GetXaxis()->SetTitle("#eta_{rec}/#eta_{gen}/");
215 Etahist[indexH]->GetYaxis()->SetTitle("#");
216 }
217
218 //FOR HF
219 EGen_hf = new TH1D( "GenJetE_hf" , "GenJet E_{T} for E_{T}>20 GeV and 3<#||{#eta}<5 ", nbinsE,Emin, Emax ) ;
220 EGen_hf->Sumw2();
221 ECalo_hf = new TH1D( "CaloJetE_hf" , "CaloJet E_{T} for QCD jets Pt_15_120 Gev & 3<#||{#eta}<5",nbinsE,Emin, Emax) ;
222 ECalo_hf->Sumw2();
223 EGen_hf_x = new TH1D( "GenJetE_hf_x" , "GenJet E_{T} for E_{T}>20 GeV and 3<#||{#eta}<5 ", nbinsE,Emin, Emax ) ;
224 EGen_hf_x->Sumw2();
225 corr_ECalo_hf = new TH1D("corr_CaloJetE_hf" , "corr_CaloJet E_{T} for E_{T}>20 GeV and 3<#||{#eta}<5",nbinsE,Emin, Emax) ;
226 corr_ECalo_hf->Sumw2();
227
228 RatioE_hf_vs_Eta= new TProfile(" RatioE_vs_Eta_hf","E^{reco}_{T}/E^{MC}_{T} for QCD jets Pt_15_120 Gev & 3<#||{#eta}<5 ", 20, 3.,5.);
229 RatioE_hf_vs_E= new TProfile("RatioE_vs_E_hf","E^{reco}_{T}/E^{MC}_{T} for QCD jets Pt_15_120 Gev & HF ", nbinsE,Emin ,Emax);
230 distX_hf = new TH2D("distX_hf","",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
231 ///
232 distX_1_hf = new TH2D("distX_1_hf"," x>0.1",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
233 distX_Matching_parton_hf = new TH2D("distX_Matchin_parton_hf"," Matching with parton cone=0.04",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
234 //FOR CASTOR
235 distX_castor = new TH2D("distX_castor","",nbinsBjoX,BjoXmin,BjoXmax,200,0. ,200.);
236 distX_Matching_parton_CASTOR = new TH2D("distX_Matchin_parton_CASTOR"," Matching with parton cone=0.04",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
237 //X dist for RecJet at CASTOR & HF
238 distX_Matching_parton_hf_RecJet = new TH2D("distX_Matchin_parton_hf_RecJet"," Matching with parton cone=0.04 ",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
239 distX_Matching_parton_CASTOR_RecJet = new TH2D("distX_Matchin_parton_CASTOR_RecJet"," Matching with parton cone=0.04",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
240 return ;
241 }
242
243
244
245 void ForwardJet::analyze( const Event& e, const EventSetup& ){
246
247 //Constant
248 double PtCut=20.;//GeV
249 double ETCut=20; //GeV
250 double threshold = 0.05;
251 double threshold_1 = 0.04;// min distance between parton and GenJEt
252 double sqrts = 14000.;//Gev
253 float BjoX;//Bjorken-x for GenJet
254 float BjoX_Rec;//Bjorken-x for RecJet
255 double limBins[] = {10,20,30,40,50,60,70,80,100.0,120.,150.,200};
256 int nbreBins=11;
257
258 float min_distance=1.;//maximum values of minumum distance between parton and GenJEt
259 //cout << "Event No: " << e.id() << " started to be analyzed" << endl;
260 if(eventnum%100==0)cout<<"EventNumber====>"<<eventnum<<endl;
261
262 ///iterativecone5
263 Handle< GenJetCollection > GenJetsHandle ;
264 try{
265 e.getByLabel( "iterativeCone5GenJets", GenJetsHandle );
266 }catch(cms::Exception& ex){
267 LogError("ForwardJet") << "Error! can't get collection.";
268 }
269
270 Handle< CaloJetCollection > CaloJetsHandle ;
271 try{
272 e.getByLabel( "iterativeCone5CaloJets", CaloJetsHandle );
273 }catch(cms::Exception& ex){
274 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
275 }
276 // Handle<CaloJetCollection> CorJetsHandle;
277 //try{
278 // e.getByLabel( CorJetAlgorithm, CorJetsHandle );
279 //}catch(cms::Exception& ex){
280 //LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
281 //}
282
283 //
284 //ktJets
285 /* Handle< GenJetCollection > GenJetsHandle ;
286 try{
287 e.getByLabel( "Kt10GenJets", GenJetsHandle );
288 }catch(cms::Exception& ex){
289 LogError("ForwardJet") << "Error! can't get collection.";
290 }
291
292 Handle< CaloJetCollection > CaloJetsHandle ;
293 try{
294 e.getByLabel( "Kt10CaloJets", CaloJetsHandle );
295 }catch(cms::Exception& ex){
296 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
297 }*/
298
299
300 ///MidPoint
301 /*Handle< GenJetCollection > GenJetsHandle ;
302 try{
303 e.getByLabel( "midPointCone5GenJets", GenJetsHandle );
304 }catch(cms::Exception& ex){
305 LogError("ForwardJet") << "Error! can't get collection.";
306 }
307
308 Handle< CaloJetCollection > CaloJetsHandle ;
309 try{
310 e.getByLabel( "midPointCone5CaloJets", CaloJetsHandle );
311 }catch(cms::Exception& ex){
312 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
313 }*/
314
315 /*************************************************/
316 // GEN JET //
317 /*************************************************/
318
319
320 vector<GenJet> tmpGenJet;
321 tmpGenJet.clear();
322 if(GenJetsHandle->size())
323 {
324 for(GenJetCollection::const_iterator it=GenJetsHandle->begin();it!=GenJetsHandle->end();it++)
325 tmpGenJet.push_back(*it);
326 std::stable_sort(tmpGenJet.begin(),tmpGenJet.end(),GenJetSort());
327 }
328
329 GenJet fGJ;
330 if(tmpGenJet.size()>0){
331 for(int g = 0; g <tmpGenJet.size();++g){
332
333 fGJ = tmpGenJet[g];
334 //fGJ.print();
335 //cout<<"particle Id=="<<fGJ.pdgId()<<endl;
336 if(fGJ.et()>ETCut)
337 {
338 ptGen->Fill(fGJ.pt());
339 yGen->Fill(fGJ.eta());
340 phiGen->Fill(fGJ.phi());
341 EGen->Fill(fGJ.et());
342 }
343 if (abs(fGJ.eta())>=3. && abs(fGJ.eta())<=5.)//Distribution BjoX and E for HF
344 {
345
346 if(fGJ.et()>ETCut)
347 {
348
349 BjoX=log10((fGJ.pt()/sqrts)*exp((-1)*(fGJ.eta())));
350 distX_hf->Fill( BjoX,fGJ.et());
351 EGen_hf->Fill(fGJ.et());
352
353 if(BjoX>=-1)
354 {
355 distX_1_hf->Fill( BjoX,fGJ.et());
356 }
357
358 }
359 }
360
361 if (abs(fGJ.eta())>=5.3 && abs(fGJ.eta())<=6.6)//Distribution BjoX for CASTOR
362 {
363
364 if(fGJ.et()>ETCut)
365 {
366 // cout<<" Castor Eta Bjr=="<<fGJ.eta()<< " GenEt=="<<fGJ.eta()<<endl;
367 BjoX=log10((fGJ.pt()/sqrts)*exp((-1)*(fGJ.eta())));
368 distX_castor->Fill( BjoX,fGJ.et());
369
370
371
372 }
373 }
374
375 }
376 }
377 /*******************************************/
378 // Matching parton dR<0.04 //
379 /******************************************/
380
381 GenJet fGenJet;
382
383 jfi.readEvent(e);//for parton
384
385 //cout <<"\nList of partons:\n";
386 vector<MCParton> partonList = jfi.getListOfPartons();
387 //Fill partonPt
388 if(partonList.size()>0){
389 for (vector<MCParton>::iterator i = partonList.begin(); i != partonList.end(); ++i)
390 {
391 float partonEta=i->fourVector().eta();
392
393 if (abs(partonEta)>=3. && abs(partonEta)<=5.)//Distribution for HF
394 {
395 if ((i->fourVector().pt())>PtCut)
396 {
397 ptParton->Fill(i->fourVector().pt());
398 }
399 }
400 }
401 }
402 //corrCaloJet
403 Handle<CaloJetCollection> corjets;
404 try{
405 e.getByLabel(CorJetAlgorithm, corjets );
406 }catch(cms::Exception& ex){
407 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
408 }
409
410 if(tmpGenJet.size()>0 && partonList.size()>0)
411 {
412
413 //cout << "Found " << partonList.size() << " partons" << endl;
414 //cout << "Found " <<tmpGenJet.size() << " genjet" << endl;
415
416
417 for(int gg = 0; gg <tmpGenJet.size();++gg)
418 {
419 fGenJet= tmpGenJet[gg];
420
421 float GenJETEta=fGenJet.eta();
422 float GenJETPhi=fGenJet.phi();
423
424
425 for (vector<MCParton>::iterator i = partonList.begin(); i != partonList.end(); ++i)
426 {
427 //cout << "Found Parton of flavour: "<< i->flavour() << endl;
428 //cout << " Eta/phi of parton: " << i->summedDaughterMomentum().eta()
429 //<<" , "<< i->summedDaughterMomentum().phi()<<endl;
430
431 float PartonEta=i->fourVector().eta();
432 float PartonPhi=i->fourVector().phi();
433 float PartonPt=i->fourVector().pt();
434
435
436 for(CaloJetCollection::const_iterator cor=corjets->begin();cor!=corjets->end();cor++)
437 {
438
439 float CaloJetEta=cor->eta();
440 float CaloJetPhi=cor->phi();
441
442
443 if (abs(fGenJet.eta())>=3. && abs(fGenJet.eta())<=5.)//Distribution for HF
444 {
445 float distance_1 = delR(GenJETEta ,PartonEta,GenJETPhi,PartonPhi);//distance for Parton&GenJet
446
447
448
449
450 //*******************************************//
451 //Find Min distance between parton and GenJet//
452 //*******************************************//
453
454 if (distance_1 < min_distance)
455 {
456 min_distance=distance_1 ;
457 }
458 Min_deltaR=min_distance;
459
460 //Fill GenJet_Et-PartonEt for Different DeltaR
461
462 if (fGenJet.pt()>PtCut)
463 {
464 if(distance_1<0.02){ Min_DeltaR_02->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
465 if(distance_1<0.04){ Min_DeltaR_04->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
466 if(distance_1<0.06){ Min_DeltaR_06->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
467 if(distance_1<0.08){ Min_DeltaR_08->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
468 if(distance_1<0.1){ Min_DeltaR_1->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
469 if(distance_1<0.2){ Min_DeltaR_2->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
470 if(distance_1<0.5){ Min_DeltaR_5->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
471 if(distance_1<0.8){ Min_DeltaR_8->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
472 }
473
474 //*******************************************//
475
476
477 if(fGenJet.et()>ETCut && distance_1<threshold_1)
478 {
479
480 BjoX=log10((fGenJet.pt()/sqrts)*exp((-1)*(fGenJet.eta())));
481
482 if(BjoX>=-1)//x>0.1
483 {
484 pid_highx_parton_01_hf->Fill(i->flavour());
485 }
486 if(BjoX>=-0.698)//X>0.2
487 {
488 pid_highx_parton_02_hf->Fill(i->flavour());
489 }
490 if(BjoX>=-0.397)//x>0.4
491 {
492 pid_highx_parton_04_hf->Fill(i->flavour());
493 }
494 if(BjoX>=-0.222)//x>0.6
495 {
496 pid_highx_parton_06_hf->Fill(i->flavour());
497 }
498
499 distX_Matching_parton_hf->Fill( BjoX,fGenJet.et());
500 pid_Allx_parton_hf->Fill(i->flavour());
501 EGen_hf_x->Fill(fGenJet.et());
502
503 }
504
505 /********************************************/
506 // Matching GenJet&Parton&RecJet for HF //
507 /*******************************************/
508 float distance_2 = delR(GenJETEta ,CaloJetEta,GenJETPhi,CaloJetPhi);//distance for GenJet&RecJet
509 if(distance_1<threshold_1 && distance_2<threshold && cor->et()>ETCut)
510 {
511
512 BjoX_Rec=log10((cor->pt()/sqrts)*exp((-1)*(cor->eta())));
513
514 distX_Matching_parton_hf_RecJet->Fill(BjoX_Rec,fGenJet.et());
515
516 }
517 }//HF eta
518
519
520 if (abs(fGenJet.eta())>=5.3 && abs(fGenJet.eta())<=6.6)//for CASTOR
521 {
522
523 float distance_1 = delR(GenJETEta ,PartonEta,GenJETPhi,PartonPhi);//distance for Parton&GenJet
524
525
526 if(fGenJet.et()>ETCut && distance_1<threshold_1)
527 {
528
529 BjoX=log10((fGenJet.pt()/sqrts)*exp((-1)*(fGenJet.eta())));
530
531 if(BjoX>=-1)//x>0.1
532 {
533 pid_highx_parton_01_castor->Fill(i->flavour());
534 }
535 if(BjoX>=-0.698)//X>0.2
536 {
537 pid_highx_parton_02_castor->Fill(i->flavour());
538 }
539 if(BjoX>=-0.397)//x>0.4
540 {
541 pid_highx_parton_04_castor->Fill(i->flavour());
542 }
543 if(BjoX>=-0.222)//x>0.6
544 {
545 pid_highx_parton_06_castor->Fill(i->flavour());
546 }
547
548 distX_Matching_parton_CASTOR->Fill( BjoX,fGenJet.et());
549 pid_Allx_parton_CASTOR->Fill(i->flavour());
550
551
552
553
554 }
555 /********************************************/
556 // Matching GenJet&Parton&RecJet for CASTOR //
557 /*******************************************/
558 float distance_2 = delR(GenJETEta ,CaloJetEta,GenJETPhi,CaloJetPhi);//distance for GenJet&RecJet
559 if(distance_1<threshold_1 && distance_2<threshold && cor->et()>ETCut)
560 {
561
562 BjoX_Rec=log10((cor->pt()/sqrts)*exp((-1)*(cor->eta())));
563
564 distX_Matching_parton_CASTOR_RecJet->Fill(BjoX_Rec,fGenJet.et());
565
566 }
567
568
569 }//CASTOR eta
570
571 }//correct RecJet
572 }//Parton
573
574 }//GenJEt
575
576 }
577
578
579 Min_DeltaR->Fill(Min_deltaR);
580
581 /*************************************************/
582 // REC0 JET //
583 /*************************************************/
584
585 vector<CaloJet> tmpCaloJet;
586 tmpCaloJet.clear();
587 if(CaloJetsHandle->size()){
588 for(CaloJetCollection::const_iterator it=CaloJetsHandle->begin();it!=CaloJetsHandle->end();it++)
589 tmpCaloJet.push_back(*it);
590 std::stable_sort(tmpCaloJet.begin(),tmpCaloJet.end(),CaloJetSort());
591 }
592 //cout << "There are " << tmpCaloJet.size() << " CaloJets" <<endl;
593
594 // CaloJet TotRecJet;
595 CaloJet fRJ;
596 if(tmpCaloJet.size()>0){
597 for(int r = 0; r <tmpCaloJet.size();++r)
598 {
599 fRJ = tmpCaloJet[r];
600 if(fRJ.et()>ETCut)
601 {
602 ptCalo->Fill(fRJ.pt());
603 yCalo->Fill(fRJ.eta());
604 phiCalo->Fill(fRJ.phi());
605 ECalo->Fill(fRJ.et());
606 }
607
608 if (abs(fRJ.eta())>=3. && abs(fRJ.eta())<=5.)
609 {
610 if(fRJ.et()>ETCut)
611 {
612
613 ECalo_hf->Fill(fRJ.et());
614 }
615 }
616
617 }
618 }
619
620 /*************************************************/
621 // CORRECTED REC0 JET //
622 /*************************************************/
623 Handle<CaloJetCollection> CorJets;
624 try{
625 e.getByLabel(CorJetAlgorithm, CorJets );
626 }catch(cms::Exception& ex){
627 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
628 }
629
630 for(CaloJetCollection::const_iterator corjets=CorJets->begin();corjets!=CorJets->end();corjets++)
631 {
632
633 if(corjets->et()>ETCut)
634 {
635 corr_ptCalo->Fill(corjets->pt());
636 corr_ECalo->Fill(corjets->et());
637 }
638
639 if (abs(corjets->eta())>=3. && abs(corjets->eta())<=5.)
640 {
641 if(corjets->et()>ETCut)
642 {
643
644 corr_ECalo_hf->Fill(corjets->et());
645 }
646 }
647 }
648
649
650
651 /*******************************************/
652 // Matching dR<0.2 //
653 /******************************************/
654
655 GenJet fGJet;
656 if(tmpGenJet.size()>0 )
657 {
658 for(int gg = 0; gg <tmpGenJet.size();++gg)
659 {
660 fGJet = tmpGenJet[gg];
661
662 float GenJetEta=fGJet.eta();
663 float GenJetPhi=fGJet.phi();
664
665
666 Handle<CaloJetCollection> corjets;
667 try{
668 e.getByLabel(CorJetAlgorithm, corjets );
669 }catch(cms::Exception& ex){
670 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
671 }
672
673 for(CaloJetCollection::const_iterator cor=corjets->begin();cor!=corjets->end();cor++)
674 {
675
676 float CaloJetEta=cor->eta();
677 float CaloJetPhi=cor->phi();
678
679 float distance = delR(GenJetEta ,CaloJetEta,GenJetPhi,CaloJetPhi);
680
681 if (distance<=threshold && fGJet.et()>ETCut )
682 {
683
684 RatioE_vs_Eta->Fill(abs(GenJetEta), (cor->et())/fGJet.et());
685 RatioE_vs_E->Fill(fGJet.et(), (cor->et())/fGJet.et());
686
687 if (abs(cor->eta())>=3. && abs(cor->eta())<=5.)// FOR HF
688 {
689 double ESim= fGJet.et();
690 int trueIndex=determineIndex(limBins,ESim,nbreBins);
691 if(trueIndex!=-1)
692 {
693
694
695 EThist[trueIndex]->Fill((cor->et())/(fGJet.et()));
696 Phihist[trueIndex]->Fill(abs(cor->phi())/abs(fGJet.phi()));
697 Etahist[trueIndex]->Fill(abs(cor->eta())/abs(fGJet.eta()));
698 trueRecoJet->Fill(cor->et());
699 countD++;
700 }
701 RatioE_hf_vs_Eta->Fill(abs(GenJetEta), (cor->et())/fGJet.et());
702 RatioE_hf_vs_E->Fill(fGJet.et(),(cor->et()) /fGJet.et());
703
704 }
705 }
706
707
708 // for not corrected calo jet
709
710 /* for(int rr = 0; rr <tmpCaloJet.size();++rr){
711 fRJet = tmpCaloJet[rr];
712
713 float CaloJetEta=fRJet.eta();
714 float CaloJetPhi=fRJet.phi();
715
716 float distance = delR(GenJetEta ,CaloJetEta,GenJetPhi,CaloJetPhi);
717
718 // if distance < threshold , fill histogramm
719 if (distance<=threshold && fGJet.et()>ETCut )
720 {
721
722 //if(genjetused[gg] == 0)
723 //{
724
725 RatioE_vs_Eta->Fill(abs(GenJetEta), fRJet.et()/fGJet.et());
726 RatioE_vs_E->Fill(fGJet.et(), fRJet.et()/fGJet.et());
727
728 if (abs(fRJet.eta())>=3. && abs(fRJet.eta())<=5.)// FOR HF
729 {
730 double ESim= fGJet.et();
731 int trueIndex=determineIndex(limBins,ESim,nbreBins);
732 if(trueIndex!=-1)
733 {
734 cout<<"RjetEta=="<<fRJet.eta()<<" GenJetEt=="<<fGJet.et()<<endl;
735
736 EThist[trueIndex]->Fill((fRJet.et())/(fGJet.et()));
737 Phihist[trueIndex]->Fill(abs(fRJet.phi())/abs(fGJet.phi()));
738 Etahist[trueIndex]->Fill(abs(fRJet.eta())/abs(fGJet.eta()));
739 trueRecoJet->Fill(fRJet.et());
740 countD++;
741 }
742 RatioE_hf_vs_Eta->Fill(abs(GenJetEta), fRJet.et()/fGJet.et());
743 RatioE_hf_vs_E->Fill(fGJet.et(), fRJet.et()/fGJet.et());
744
745 }
746 }
747 }*/
748 }
749 }
750 }
751
752
753 //RecGenTree->Fill();
754 //cout << endl;
755 //cout << "===> " << e.id() << " finished" << endl;
756 eventnum++;
757 return ;
758
759 }
760
761 void ForwardJet::endJob(){
762 long int sigma=146e7;//pb for pt_15_20 from crossection_QCD_event-list.pdf
763
764
765 ptGen->Scale(sigma/(eventnum*ptGen->GetBinWidth(1)));
766 ptParton->Scale(sigma/(eventnum*ptParton->GetBinWidth(1)));
767 yGen->Scale(sigma/(eventnum*yGen->GetBinWidth(1)));
768 phiGen->Scale(sigma/(eventnum*phiGen->GetBinWidth(1)));
769 EGen->Scale(sigma/(eventnum* EGen->GetBinWidth(1)));
770 ptCalo->Scale(sigma/(eventnum*ptCalo->GetBinWidth(1)));
771 corr_ptCalo->Scale(sigma/(eventnum*corr_ptCalo->GetBinWidth(1)));
772 yCalo->Scale(sigma/(eventnum*yCalo->GetBinWidth(1)));
773 phiCalo->Scale(sigma/(eventnum*phiCalo->GetBinWidth(1)));
774 ECalo->Scale(sigma/(eventnum*ECalo->GetBinWidth(1)));
775 corr_ECalo->Scale(sigma/(eventnum*corr_ECalo->GetBinWidth(1)));
776 EGen_hf->Scale(sigma/(eventnum*EGen_hf->GetBinWidth(1)));
777 EGen_hf_x->Scale(sigma/(eventnum*EGen_hf_x->GetBinWidth(1)));
778 ECalo_hf->Scale(sigma/(eventnum*ECalo_hf->GetBinWidth(1)));
779 corr_ECalo_hf->Scale(sigma/(eventnum*corr_ECalo_hf->GetBinWidth(1)));
780
781 fOutputFile->Write() ;
782 fOutputFile->Close() ;
783 return ;
784 }
785 float ForwardJet::delR(const float& eta1,const float& eta2,const float& phi1,const float& phi2)
786 {
787
788 float Distance =sqrt((2*atan(tan((phi1-phi2)/2)))*(2*atan(tan((phi1-phi2)/2)))+(eta1-eta2)*(eta1-eta2));
789 return Distance;
790
791 }
792 int ForwardJet::determineIndex( double* limBins, double data, int nbreBins){
793
794 int result =-1;
795 if(data<limBins[0]) return result;
796
797 for(int i = 0; i<nbreBins;i++){
798 if(limBins[i]<= data && data<limBins[i+1]){
799 result=i;
800 break;}
801 result=i;
802 }
803
804 return result;
805
806 }
807
808 DEFINE_FWK_MODULE(ForwardJet);