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

# User Rev Content
1 salim 1.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 salim 1.2 ForwardJet::ForwardJet( const ParameterSet& pset ) : fOutputFileName( pset.getUntrackedParameter<string>("HistOutFile",string("/data2/salim/ForwardJet/ForwardJet_QCD_pt_120_170.root"))),fOutputFile(0)
63 salim 1.1
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 salim 1.2 int nbinspt = 300;
77 salim 1.1 int nbinsphi = 100;
78 salim 1.2 int nbinsE= 300;
79 salim 1.1 int nbinsEta =100;
80     int nbinsBjoX=600;
81     double ptmin = 0; //GeV
82 salim 1.2 double ptmax = 300; //GeV
83 salim 1.1 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 salim 1.2 double Emax=300; //GeV
91 salim 1.1 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 salim 1.2 double limBins[] = {10,20,30,40,50,60,70,80,100.0,120.,150.,250};
133 salim 1.1
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 salim 1.2 sprintf(nameHist1,"phidistribution%i",indexH);
147     sprintf(nameHist2,"etadistribution%i",indexH);
148 salim 1.1
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 salim 1.2 Phihist[indexH]=new TH1D(nameHist1,label1,480, 0.2, 1.8 );
159 salim 1.1 Phihist[indexH]->Sumw2();
160 salim 1.2 Phihist[indexH]->GetXaxis()->SetTitle("#phi_{rec}/#phi_{gen}");
161 salim 1.1 Phihist[indexH]->GetYaxis()->SetTitle("#");
162    
163 salim 1.2 Etahist[indexH]=new TH1D(nameHist2,label2,400,0.8, 1.2);
164 salim 1.1 Etahist[indexH]->Sumw2();
165 salim 1.2 Etahist[indexH]->GetXaxis()->SetTitle("#eta_{rec}/#eta_{gen}/");
166 salim 1.1 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 salim 1.2 //cout << "Event No: " << e.id() << " started to be analyzed" << endl;
198 salim 1.1 //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 salim 1.2 //fGJ.print();
268     //cout<<"particle Id=="<<fGJ.pdgId()<<endl;
269 salim 1.1 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 salim 1.2 // cout<<" Castor Eta Bjr=="<<fGJ.eta()<< " GenEt=="<<fGJ.eta()<<endl;
304 salim 1.1 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 salim 1.2 // for(int gg = 0; gg <500;++gg){genjetused[gg] = 0;}
363 salim 1.1
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 salim 1.2 //{
388 salim 1.1
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 salim 1.2 Phihist[trueIndex]->Fill(abs(fRJet.phi())/abs(fGJet.phi()));
402     Etahist[trueIndex]->Fill(abs(fRJet.eta())/abs(fGJet.eta()));
403 salim 1.1 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 salim 1.2 //}
410     }
411 salim 1.1
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);