ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SCerci/src/ForwardJet.cc
Revision: 1.2
Committed: Wed Jul 25 19:40:27 2007 UTC (17 years, 9 months ago) by salim
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +22 -22 lines
Log Message:
*** empty log message ***

File Contents

# Content
1 /////////////////////////////////////////////////////
2 // //
3 // Original Author: Salim Cerci //
4 // Created: Mon Jul 16 18:50:23 CET 2007 //
5 // //
6 // //
7 /////////////////////////////////////////////////////
8 #include <iostream>
9 #include "Analysis/ForwardJet/interface/ForwardJet.h"
10
11 #include "DataFormats/Common/interface/Ref.h"
12 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
13 #include "DataFormats/Common/interface/Ref.h"
14 #include "DataFormats/JetReco/interface/Jet.h"
15 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
16 #include "DataFormats/JetReco/interface/BasicJetCollection.h"
17 #include "DataFormats/JetReco/interface/GenJet.h"
18 #include "DataFormats/JetReco/interface/CaloJet.h"
19 #include "DataFormats/METReco/interface/CaloMET.h"
20 #include "DataFormats/METReco/interface/GenMET.h"
21 #include "DataFormats/JetReco/interface/GenJetfwd.h"
22 #include "DataFormats/JetReco/interface/BasicJet.h"
23 #include "DataFormats/JetReco/interface/BasicJetfwd.h"
24 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
25
26 #include "CLHEP/HepMC/ReadHepMC.h"
27 #include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
28
29 #include "TFile.h"
30 #include "TTree.h"
31 #include "TH1.h"
32 #include "TH1F.h"
33 #include "TH2F.h"
34 #include "TMath.h"
35 #include "TF1.h"
36 #include "TCanvas.h"
37 #include "TGraphErrors.h"
38 #include "TProfile.h"
39
40 using namespace edm;
41 using namespace std;
42 using namespace reco;
43
44
45 int eventnum=1;
46 class GenJetSort{
47 public:
48 bool operator()(const GenJet& a, const GenJet& b) {
49 return a.et() > b.et();
50 }
51 };
52
53 class CaloJetSort{
54 public:
55 bool operator()(const CaloJet& a, const CaloJet& b) {
56 return a.et() > b.et();
57 }
58 };
59 //ForwardJet::ForwardJet( const ParameterSet& pset )
60 //{}
61
62 ForwardJet::ForwardJet( const ParameterSet& pset ) : fOutputFileName( pset.getUntrackedParameter<string>("HistOutFile",string("/data2/salim/ForwardJet/ForwardJet_QCD_pt_120_170.root"))),fOutputFile(0)
63
64 {
65 }
66
67 void ForwardJet::beginJob( const EventSetup& )
68 {
69 fOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
70
71 //RecGenTree = new TTree("RecJet_GenJet", "Analysis");
72 //RecGenTree->Branch("fNRecJ", &fNRecJ, "fNRecJ/I");
73 //RecGenTree->Branch("fRecJetPt", &fRecJetPt, "fRecJetPt[fNRecJ]/D");
74
75 // limits for histograms
76 int nbinspt = 300;
77 int nbinsphi = 100;
78 int nbinsE= 300;
79 int nbinsEta =100;
80 int nbinsBjoX=600;
81 double ptmin = 0; //GeV
82 double ptmax = 300; //GeV
83 double ymin = -10;
84 double ymax = 10;
85 double BjoXmin = -6;
86 double BjoXmax = 0;
87 double phimin = -TMath::Pi();
88 double phimax = TMath::Pi();
89 double Emin=0; //GeV
90 double Emax=300; //GeV
91 double ETCut=20; //GeV
92 double threshold = 0.2;
93 double sqrts = 14000.;//Gev
94 float BjoX;//Bjorken-x
95 // create histograms
96 ptGen = new TH1D("ptGen","Gen jets pT for E_{T}>20 GeV ",nbinspt,ptmin,ptmax);
97 ptGen->Sumw2();
98 yGen = new TH1D("yGen","Gen jets rapidity for E_{T}>20 GeV ",nbinsEta,ymin,ymax);
99 yGen->Sumw2();
100 phiGen = new TH1D("phiGen","Gen jets phi for E_{T}>20 GeV ",nbinsphi,phimin,phimax);
101 phiGen->Sumw2();
102 EGen = new TH1D ("EGen","Gen jets energy for E_{T}>20 GeV ",nbinsE,Emin, Emax);
103 EGen->Sumw2();
104
105 ptCalo = new TH1D("ptCalo","Calo jets pT for E_{T}>20 GeV ",nbinspt,ptmin,ptmax);
106 ptCalo->Sumw2();
107 yCalo = new TH1D("yCalo","Calo jets rapidity for E_{T}>20 GeV ",nbinsEta,ymin,ymax);
108 yCalo->Sumw2();
109 phiCalo = new TH1D("phiCalo","Calo jets phi for E_{T}>20 GeV ",nbinsphi,phimin,phimax);
110 phiCalo->Sumw2();
111 ECalo = new TH1D ("ECalo","Calo jets energy for E_{T}>20 GeV ",nbinsE,Emin, Emax);
112 ECalo->Sumw2();
113
114 trueRecoJet = new TH1D ("trueRecoJet","True Reco Jet for E_{T}>20 GeV ",nbinsE,Emin, Emax);
115 trueRecoJet->Sumw2();
116
117 RatioE_vs_Eta= new TProfile(" RatioE_vs_Eta","E^{reco}_{T}/E^{MC}_{T} for E_{T}>20 GeV ", (nbinsEta/2), 0.,5.);
118 RatioE_vs_E= new TProfile("RatioE_vs_E ","E^{reco}_{T}/E^{MC}_{T} for E_{T}>20>20 GeV ", nbinsE,Emin ,Emax);
119
120 pid_highx_parton = new TH1D("pid_highx_parton","PID of parton (jet)with x>0.1 ",30,-0.5,29.5);
121
122 //Creating a binning in ET and an histogram for each of the ET bin
123 unsigned int nbreBins=11;
124 //EThist [11];
125 char nameHist [50];//for Et distribution
126 char nameHist1 [50];//for reco jet phi distribution
127 char nameHist2 [50];//for reco jet eta distribution
128 char label [500];//for Et distribution
129 char label1 [500];//for reco jet phi distribution
130 char label2 [500];//for reco jet eta distribution
131
132 double limBins[] = {10,20,30,40,50,60,70,80,100.0,120.,150.,250};
133
134 ratioE = new TH1D ("ratioE", "E^{reco}_{T}/E^{MC}_{T}",nbreBins, Emin, Emax);
135 ratioE->Sumw2();
136 ratioE->SetBins(nbreBins,limBins);
137
138 sigmaErEm = new TH1D ("sigmaErEm", "sigma Ereco / Emc",nbreBins, Emin, Emax);
139
140 sigmaErEm->Sumw2();
141 sigmaErEm->SetBins(nbreBins,limBins);
142
143 for (unsigned int indexH=0;indexH<nbreBins;indexH++)
144 {
145 sprintf(nameHist,"ETdistribution%i",indexH);
146 sprintf(nameHist1,"phidistribution%i",indexH);
147 sprintf(nameHist2,"etadistribution%i",indexH);
148
149 sprintf(label,"E_{T}(reco)/E_{T}(Gen) #in #left[%f , %f #right] GeV",limBins[indexH],limBins[indexH+1]);
150 sprintf(label1,"#phi_{rec} #in #left[%f , %f #right] GeV",limBins[indexH],limBins[indexH+1]);
151 sprintf(label2,"#eta_{rec} #in #left[%f , %f #right] GeV",limBins[indexH],limBins[indexH+1]);
152
153 EThist[indexH]=new TH1D(nameHist,label,nbinsE,0., 1.3);
154 EThist[indexH]->Sumw2();
155 EThist[indexH]->GetXaxis()->SetTitle("E_{T}^{reco}/E_{t}^{gen}");
156 EThist[indexH]->GetYaxis()->SetTitle("#");
157
158 Phihist[indexH]=new TH1D(nameHist1,label1,480, 0.2, 1.8 );
159 Phihist[indexH]->Sumw2();
160 Phihist[indexH]->GetXaxis()->SetTitle("#phi_{rec}/#phi_{gen}");
161 Phihist[indexH]->GetYaxis()->SetTitle("#");
162
163 Etahist[indexH]=new TH1D(nameHist2,label2,400,0.8, 1.2);
164 Etahist[indexH]->Sumw2();
165 Etahist[indexH]->GetXaxis()->SetTitle("#eta_{rec}/#eta_{gen}/");
166 Etahist[indexH]->GetYaxis()->SetTitle("#");
167 }
168
169 //FOR HF
170 EGen_hf = new TH1D( "GenJetE_hf" , "GenJet E_{T} for E_{T}>20 GeV and 3<#||{#eta}<5 ", nbinsE,Emin, Emax ) ;
171 EGen_hf->Sumw2();
172 ECalo_hf = new TH1D( "CaloJetE_hf" , "CaloJet E_{T} for E_{T}>20 GeV and 3<#||{#eta}<5",nbinsE,Emin, Emax) ;
173 ECalo_hf->Sumw2();
174 RatioE_hf_vs_Eta= new TProfile(" RatioE_vs_Eta_hf","E^{reco}_{T}/E^{MC}_{T} for E_{T}>20 GeV and 3<#||{#eta}<5 ", 20, 3.,5.);
175 RatioE_hf_vs_E= new TProfile("RatioE_vs_E_hf","E^{reco}_{T}/E^{MC}_{T} for E_{T}>20 GeV and hf ", nbinsE,Emin ,Emax);
176 distX_hf = new TH2D("distX_hf","",nbinsBjoX,BjoXmin,BjoXmax,nbinsE,Emin ,Emax);
177 ///
178 checkdistX_hf = new TH1D("distX_chech_hf","",100,-4.,-2.);
179 //FOR CASTOR
180 distX_castor = new TH2D("distX_castor","",nbinsBjoX,BjoXmin,BjoXmax,nbinsE,Emin ,Emax);
181
182 return ;
183 }
184
185
186
187 void ForwardJet::analyze( const Event& e, const EventSetup& ){
188
189 //Constant
190 double ETCut=20; //GeV
191 double threshold = 0.2;
192 double sqrts = 14000.;//Gev
193 float BjoX;//Bjorken-x
194 double limBins[] = {10,20,30,40,50,60,70,80,100.0,120.,150.,200};
195 int nbreBins=11;
196
197 //cout << "Event No: " << e.id() << " started to be analyzed" << endl;
198 //cout<<"EventNumber====>"<<e.entries()<<endl;
199
200 cout<<"Event Number=="<<eventnum<<endl;
201 ///iterativecone5
202 Handle< GenJetCollection > GenJetsHandle ;
203 try{
204 e.getByLabel( "iterativeCone5GenJets", GenJetsHandle );
205 }catch(cms::Exception& ex){
206 LogError("ForwardJet") << "Error! can't get collection.";
207 }
208
209 Handle< CaloJetCollection > CaloJetsHandle ;
210 try{
211 e.getByLabel( "iterativeCone5CaloJets", CaloJetsHandle );
212 }catch(cms::Exception& ex){
213 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
214 }
215
216 //
217 //ktJets
218 /* Handle< GenJetCollection > GenJetsHandle ;
219 try{
220 e.getByLabel( "Kt10GenJets", GenJetsHandle );
221 }catch(cms::Exception& ex){
222 LogError("ForwardJet") << "Error! can't get collection.";
223 }
224
225 Handle< CaloJetCollection > CaloJetsHandle ;
226 try{
227 e.getByLabel( "Kt10CaloJets", CaloJetsHandle );
228 }catch(cms::Exception& ex){
229 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
230 }*/
231
232
233 ///MidPoint
234 /*Handle< GenJetCollection > GenJetsHandle ;
235 try{
236 e.getByLabel( "midPointCone5GenJets", GenJetsHandle );
237 }catch(cms::Exception& ex){
238 LogError("ForwardJet") << "Error! can't get collection.";
239 }
240
241 Handle< CaloJetCollection > CaloJetsHandle ;
242 try{
243 e.getByLabel( "midPointCone5CaloJets", CaloJetsHandle );
244 }catch(cms::Exception& ex){
245 LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
246 }*/
247
248 /*************************************************/
249 // GEN JET //
250 /*************************************************/
251
252
253 vector<GenJet> tmpGenJet;
254 tmpGenJet.clear();
255 if(GenJetsHandle->size())
256 {
257 for(GenJetCollection::const_iterator it=GenJetsHandle->begin();it!=GenJetsHandle->end();it++)
258 tmpGenJet.push_back(*it);
259 std::stable_sort(tmpGenJet.begin(),tmpGenJet.end(),GenJetSort());
260 }
261
262 GenJet fGJ;
263 if(tmpGenJet.size()>0){
264 for(int g = 0; g <tmpGenJet.size();++g){
265
266 fGJ = tmpGenJet[g];
267 //fGJ.print();
268 //cout<<"particle Id=="<<fGJ.pdgId()<<endl;
269 if(fGJ.et()>ETCut)
270 {
271 ptGen->Fill(fGJ.pt());
272 yGen->Fill(fGJ.eta());
273 phiGen->Fill(fGJ.phi());
274 EGen->Fill(fGJ.et());
275 }
276 if (abs(fGJ.eta())>=3. && abs(fGJ.eta())<=5.)//Distribution BjoX and E for HF
277 {
278
279 if(fGJ.et()>ETCut)
280 {
281
282
283 BjoX=log10((fGJ.pt()/sqrts)*exp((-1)*(fGJ.eta())));
284
285
286 if(BjoX>=-1)
287 {
288 // cout<<"particle Id=="<<fGJ.pdgId()<<endl;
289 pid_highx_parton->Fill(fGJ.pdgId()); //
290 }
291
292 distX_hf->Fill( BjoX,fGJ.et());
293 EGen_hf->Fill(fGJ.et());
294
295 }
296 }
297
298 if (abs(fGJ.eta())>=5.3 && abs(fGJ.eta())<=6.6)//Distribution BjoX for CASTOR
299 {
300
301 if(fGJ.et()>ETCut)
302 {
303 // cout<<" Castor Eta Bjr=="<<fGJ.eta()<< " GenEt=="<<fGJ.eta()<<endl;
304 BjoX=log10((fGJ.pt()/sqrts)*exp((-1)*(fGJ.eta())));
305 distX_castor->Fill( BjoX,fGJ.et());
306
307 }
308 }
309
310 }
311 }
312
313
314
315 /*************************************************/
316 // REC0 JET //
317 /*************************************************/
318
319 vector<CaloJet> tmpCaloJet;
320 tmpCaloJet.clear();
321 if(CaloJetsHandle->size()){
322 for(CaloJetCollection::const_iterator it=CaloJetsHandle->begin();it!=CaloJetsHandle->end();it++)
323 tmpCaloJet.push_back(*it);
324 std::stable_sort(tmpCaloJet.begin(),tmpCaloJet.end(),CaloJetSort());
325 }
326 //cout << "There are " << tmpCaloJet.size() << " CaloJets" <<endl;
327
328 // CaloJet TotRecJet;
329 CaloJet fRJ;
330 if(tmpCaloJet.size()>0){
331 for(int r = 0; r <tmpCaloJet.size();++r)
332 {
333 fRJ = tmpCaloJet[r];
334 if(fRJ.et()>ETCut)
335 {
336 ptCalo->Fill(fRJ.pt());
337 yCalo->Fill(fRJ.eta());
338 phiCalo->Fill(fRJ.phi());
339 ECalo->Fill(fRJ.et());
340 }
341
342 if (abs(fRJ.eta())>=3. && abs(fRJ.eta())<=5.)
343 {
344 if(fRJ.et()>ETCut)
345 {
346
347 ECalo_hf->Fill(fRJ.et());
348 }
349 }
350
351 }
352 }
353
354
355 /*******************************************/
356 // Matching dR<0.2 //
357 /******************************************/
358
359 GenJet fGJet;
360 CaloJet fRJet;
361 //int genjetused[500];
362 // for(int gg = 0; gg <500;++gg){genjetused[gg] = 0;}
363
364 if(tmpGenJet.size()>0 && tmpCaloJet.size()>0)
365 {
366 for(int gg = 0; gg <tmpGenJet.size();++gg)
367 {
368 fGJet = tmpGenJet[gg];
369
370 float GenJetEta=fGJet.eta();
371 float GenJetPhi=fGJet.phi();
372
373
374 for(int rr = 0; rr <tmpCaloJet.size();++rr){
375 fRJet = tmpCaloJet[rr];
376
377 float CaloJetEta=fRJet.eta();
378 float CaloJetPhi=fRJet.phi();
379
380 float distance = delR(GenJetEta ,CaloJetEta,GenJetPhi,CaloJetPhi);
381
382 // if distance < threshold , fill histogramm
383 if (distance<=threshold && fGJet.et()>ETCut )
384 {
385
386 //if(genjetused[gg] == 0)
387 //{
388
389 RatioE_vs_Eta->Fill(abs(GenJetEta), fRJet.et()/fGJet.et());
390 RatioE_vs_E->Fill(fGJet.et(), fRJet.et()/fGJet.et());
391
392 if (abs(fRJet.eta())>=3. && abs(fRJet.eta())<=5.)// FOR HF
393 {
394 double ESim= fGJet.et();
395 int trueIndex=determineIndex(limBins,ESim,nbreBins);
396 if(trueIndex!=-1)
397 {
398 cout<<"RjetEta=="<<fRJet.eta()<<" GenJetEt=="<<fGJet.et()<<endl;
399
400 EThist[trueIndex]->Fill((fRJet.et())/(fGJet.et()));
401 Phihist[trueIndex]->Fill(abs(fRJet.phi())/abs(fGJet.phi()));
402 Etahist[trueIndex]->Fill(abs(fRJet.eta())/abs(fGJet.eta()));
403 trueRecoJet->Fill(fRJet.et());
404 countD++;
405 }
406 RatioE_hf_vs_Eta->Fill(abs(GenJetEta), fRJet.et()/fGJet.et());
407 RatioE_hf_vs_E->Fill(fGJet.et(), fRJet.et()/fGJet.et());
408
409 //}
410 }
411
412
413 }
414
415
416 }
417 }
418 }
419
420 //RecGenTree->Fill();
421 cout << endl;
422 cout << "===> " << e.id() << " finished" << endl;
423 eventnum++;
424 return ;
425
426 }
427
428 void ForwardJet::endJob(){
429 fOutputFile->Write() ;
430 fOutputFile->Close() ;
431 return ;
432 }
433 float ForwardJet::delR(const float& eta1,const float& eta2,const float& phi1,const float& phi2)
434 {
435
436 float Distance =sqrt((2*atan(tan((phi1-phi2)/2)))*(2*atan(tan((phi1-phi2)/2)))+(eta1-eta2)*(eta1-eta2));
437 return Distance;
438
439 }
440 int ForwardJet::determineIndex( double* limBins, double data, int nbreBins){
441
442 int result =-1;
443 if(data<limBins[0]) return result;
444
445 for(int i = 0; i<nbreBins;i++){
446 if(limBins[i]<= data && data<limBins[i+1]){
447 result=i;
448 break;}
449 result=i;
450 }
451
452 return result;
453
454 }
455
456 DEFINE_FWK_MODULE(ForwardJet);