ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SCerci/src/ForwardJet_v2.cc
Revision: 1.1
Committed: Thu Sep 27 13:29:57 2007 UTC (17 years, 7 months ago) by salim
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message:
*** empty log message ***

File Contents

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