ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/SCerci/src/ForwardJet_FastJ.cc
Revision: 1.2
Committed: Fri Oct 5 08:00:38 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
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     fOutputFile(0),
76    
77     CorJetAlgorithm( pset.getParameter<string>( "CorJetAlgorithm" ) )
78    
79    
80     {
81     jfi = JetFlavourIdentifier(pset.getParameter<edm::ParameterSet>("jetIdParameters"));
82     }
83    
84     void ForwardJet::beginJob( const EventSetup& )
85     {
86     fOutputFile = new TFile( fOutputFileName.c_str(), "RECREATE" ) ;
87    
88     //RecGenTree = new TTree("RecJet_GenJet", "Analysis");
89     //RecGenTree->Branch("fNRecJ", &fNRecJ, "fNRecJ/I");
90     //RecGenTree->Branch("fRecJetPt", &fRecJetPt, "fRecJetPt[fNRecJ]/D");
91    
92     // limits for histograms
93     int nbinspt = 200;
94     int nbinsphi = 100;
95     int nbinsE= 200;
96     int nbinsEta =100;
97     int nbinsBjoX=600;
98     double ptmin = 0; //GeV
99     double ptmax = 200; //GeV
100     double ymin = -10;
101     double ymax = 10;
102     double BjoXmin = -6;
103     double BjoXmax = 0;
104     double phimin = -TMath::Pi();
105     double phimax = TMath::Pi();
106     double Emin=0; //GeV
107     double Emax=200; //GeV
108     double ETCut=20; //GeV
109     double threshold = 0.2;
110     double threshold_1 = 0.04;// min distance between parton and GenJEt
111     double sqrts = 14000.;//Gev
112     float BjoX;//Bjorken-x(e^-eta)
113     float BjoX_;//Bjorken-x(e^eta)
114    
115     float Min_deltaR;// minumum distance between parton and GenJEt
116     // create histograms
117     ptGen = new TH1D("ptGen","QCD GenJet pT Pt_15_120 Gev",nbinspt,ptmin,ptmax);
118     ptGen->Sumw2();
119     yGen = new TH1D("yGen","QCD GenJet rapidity Pt_15_120 Gev ",nbinsEta,ymin,ymax);
120     yGen->Sumw2();
121     phiGen = new TH1D("phiGen","QCD GenJet Phi Pt_15_120 Gev",nbinsphi,phimin,phimax);
122     phiGen->Sumw2();
123     EGen = new TH1D ("EGen","QCD GenJet E Pt_15_120 Gev ",nbinsE,Emin, Emax);
124     EGen->Sumw2();
125     ptParton = new TH1D("ptPArton","QCD Parton pT ",nbinspt,ptmin,ptmax);
126     ptParton->Sumw2();
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    
224     corr_ECalo_hf = new TH1D("corr_CaloJetE_hf" , "corr_CaloJet E_{T} for E_{T}>20 GeV and 3<#||{#eta}<5",nbinsE,Emin, Emax) ;
225     corr_ECalo_hf->Sumw2();
226    
227     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.);
228     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);
229     distX_hf = new TH2D("distX_hf","",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
230     ///
231     distX_1_hf = new TH2D("distX_1_hf"," x>0.1",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
232     distX_Matching_parton_hf = new TH2D("distX_Matchin_parton_hf"," Matching with parton cone=0.04",nbinsBjoX,BjoXmin,BjoXmax,200,0.,200.);
233     correlationBjX=new TH2D("correlationBjX"," Correlation Bjorken-x",nbinsBjoX,BjoXmin,BjoXmax,nbinsBjoX,BjoXmin,BjoXmax);
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.2;
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     //FastJets
299     Handle< GenJetCollection > GenJetsHandle ;
300     try{
301     e.getByLabel( "Fastjet10GenJets", GenJetsHandle );
302     }catch(cms::Exception& ex){
303     LogError("ForwardJet") << "Error! can't get collection.";
304     }
305    
306     Handle< CaloJetCollection > CaloJetsHandle ;
307     try{
308     e.getByLabel( "Fastjet10CaloJets", CaloJetsHandle );
309     }catch(cms::Exception& ex){
310     LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
311     }
312    
313    
314     ///MidPoint
315     /*Handle< GenJetCollection > GenJetsHandle ;
316     try{
317     e.getByLabel( "midPointCone5GenJets", GenJetsHandle );
318     }catch(cms::Exception& ex){
319     LogError("ForwardJet") << "Error! can't get collection.";
320     }
321    
322     Handle< CaloJetCollection > CaloJetsHandle ;
323     try{
324     e.getByLabel( "midPointCone5CaloJets", CaloJetsHandle );
325     }catch(cms::Exception& ex){
326     LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
327     }*/
328    
329     /*************************************************/
330     // GEN JET //
331     /*************************************************/
332    
333    
334     vector<GenJet> tmpGenJet;
335     tmpGenJet.clear();
336     if(GenJetsHandle->size())
337     {
338     for(GenJetCollection::const_iterator it=GenJetsHandle->begin();it!=GenJetsHandle->end();it++)
339     tmpGenJet.push_back(*it);
340     std::stable_sort(tmpGenJet.begin(),tmpGenJet.end(),GenJetSort());
341     }
342    
343     GenJet fGJ;
344     if(tmpGenJet.size()>0){
345     for(int g = 0; g <tmpGenJet.size();++g){
346    
347     fGJ = tmpGenJet[g];
348     //fGJ.print();
349     //cout<<"particle Id=="<<fGJ.pdgId()<<endl;
350     if(fGJ.et()>ETCut)
351     {
352     ptGen->Fill(fGJ.pt());
353     yGen->Fill(fGJ.eta());
354     phiGen->Fill(fGJ.phi());
355     EGen->Fill(fGJ.et());
356     }
357     if (abs(fGJ.eta())>=3. && abs(fGJ.eta())<=5.)//Distribution BjoX and E for HF
358     {
359    
360     if(fGJ.et()>ETCut)
361     {
362    
363     BjoX=log10((fGJ.pt()/sqrts)*exp((-1)*(fGJ.eta())));
364     distX_hf->Fill( BjoX,fGJ.et());
365     EGen_hf->Fill(fGJ.et());
366    
367     if(BjoX>=-1)
368     {
369     distX_1_hf->Fill( BjoX,fGJ.et());
370     }
371    
372     }
373     }
374    
375     if (abs(fGJ.eta())>=5.3 && abs(fGJ.eta())<=6.6)//Distribution BjoX for CASTOR
376     {
377    
378     if(fGJ.et()>ETCut)
379     {
380     // cout<<" Castor Eta Bjr=="<<fGJ.eta()<< " GenEt=="<<fGJ.eta()<<endl;
381     BjoX=log10((fGJ.pt()/sqrts)*exp((-1)*(fGJ.eta())));
382     distX_castor->Fill( BjoX,fGJ.et());
383    
384    
385    
386     }
387     }
388    
389     }
390     }
391     /*******************************************/
392     // Matching parton dR<0.04 //
393     /******************************************/
394    
395     GenJet fGenJet;
396    
397     jfi.readEvent(e);//for parton
398    
399     //cout <<"\nList of partons:\n";
400     vector<MCParton> partonList = jfi.getListOfPartons();
401     //Fill partonPt
402     if(partonList.size()>0){
403     for (vector<MCParton>::iterator i = partonList.begin(); i != partonList.end(); ++i)
404     {
405     if ((i->fourVector().pt())>PtCut)
406     {
407     ptParton->Fill(i->fourVector().pt());
408     }
409     }
410     }
411     //corrCaloJet {
412     Handle<CaloJetCollection> corjets;
413     try{
414     e.getByLabel(CorJetAlgorithm, corjets );
415     }catch(cms::Exception& ex){
416     LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
417     }
418    
419     if(tmpGenJet.size()>0 && partonList.size()>0)
420     {
421    
422     //cout << "Found " << partonList.size() << " partons" << endl;
423     //cout << "Found " <<tmpGenJet.size() << " genjet" << endl;
424    
425    
426     for(int gg = 0; gg <tmpGenJet.size();++gg)
427     {
428     fGenJet= tmpGenJet[gg];
429    
430     float GenJETEta=fGenJet.eta();
431     float GenJETPhi=fGenJet.phi();
432    
433    
434     for (vector<MCParton>::iterator i = partonList.begin(); i != partonList.end(); ++i)
435     {
436     //cout << "Found Parton of flavour: "<< i->flavour() << endl;
437     //cout << " Eta/phi of parton: " << i->summedDaughterMomentum().eta()
438     //<<" , "<< i->summedDaughterMomentum().phi()<<endl;
439    
440     float PartonEta=i->fourVector().eta();
441     float PartonPhi=i->fourVector().phi();
442     float PartonPt=i->fourVector().pt();
443    
444    
445     for(CaloJetCollection::const_iterator cor=corjets->begin();cor!=corjets->end();cor++)
446     {
447    
448     float CaloJetEta=cor->eta();
449     float CaloJetPhi=cor->phi();
450    
451    
452     if (abs(fGenJet.eta())>=3. && abs(fGenJet.eta())<=5.)//Distribution for HF
453     {
454     float distance_1 = delR(GenJETEta ,PartonEta,GenJETPhi,PartonPhi);//distance for Parton&GenJet
455    
456    
457    
458    
459     //*******************************************//
460     //Find Min distance between parton and GenJet//
461     //*******************************************//
462    
463     if (distance_1 < min_distance)
464     {
465     min_distance=distance_1 ;
466     }
467     Min_deltaR=min_distance;
468    
469     //Fill GenJet_Et-PartonEt for Different DeltaR
470    
471     if (fGenJet.pt()>PtCut)
472     {
473     if(distance_1<0.02){ Min_DeltaR_02->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
474     if(distance_1<0.04){ Min_DeltaR_04->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
475     if(distance_1<0.06){ Min_DeltaR_06->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
476     if(distance_1<0.08){ Min_DeltaR_08->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
477     if(distance_1<0.1){ Min_DeltaR_1->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
478     if(distance_1<0.2){ Min_DeltaR_2->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
479     if(distance_1<0.5){ Min_DeltaR_5->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
480     if(distance_1<0.8){ Min_DeltaR_8->Fill(fGenJet.pt(),(fGenJet.pt()-PartonPt));}
481     }
482    
483     //*******************************************//
484    
485    
486     if(fGenJet.et()>ETCut && distance_1<threshold_1)
487     {
488    
489     BjoX=log10((fGenJet.pt()/sqrts)*exp((-1)*(fGenJet.eta())));
490     BjoX_=log10((fGenJet.pt()/sqrts)*exp((1)*(fGenJet.eta())));
491    
492     if(BjoX>=-1)//x>0.1
493     {
494     pid_highx_parton_01_hf->Fill(i->flavour());
495     }
496     if(BjoX>=-0.698)//X>0.2
497     {
498     pid_highx_parton_02_hf->Fill(i->flavour());
499     }
500     if(BjoX>=-0.397)//x>0.4
501     {
502     pid_highx_parton_04_hf->Fill(i->flavour());
503     }
504     if(BjoX>=-0.222)//x>0.6
505     {
506     pid_highx_parton_06_hf->Fill(i->flavour());
507     }
508    
509     distX_Matching_parton_hf->Fill( BjoX,fGenJet.et());
510     correlationBjX->Fill( BjoX, BjoX_);
511     pid_Allx_parton_hf->Fill(i->flavour());
512    
513     }
514    
515     /********************************************/
516     // Matching GenJet&Parton&RecJet for HF //
517     /*******************************************/
518     float distance_2 = delR(GenJETEta ,CaloJetEta,GenJETPhi,CaloJetPhi);//distance for GenJet&RecJet
519     if(distance_1<threshold_1 && distance_2<threshold && cor->et()>ETCut)
520     {
521    
522     BjoX_Rec=log10((cor->pt()/sqrts)*exp((-1)*(cor->eta())));
523    
524     distX_Matching_parton_hf_RecJet->Fill(BjoX_Rec,fGenJet.et());
525    
526     }
527     }//HF eta
528    
529    
530     if (abs(fGenJet.eta())>=5.3 && abs(fGenJet.eta())<=6.6)//for CASTOR
531     {
532    
533     float distance_1 = delR(GenJETEta ,PartonEta,GenJETPhi,PartonPhi);//distance for Parton&GenJet
534    
535    
536     if(fGenJet.et()>ETCut && distance_1<threshold_1)
537     {
538    
539     BjoX=log10((fGenJet.pt()/sqrts)*exp((-1)*(fGenJet.eta())));
540    
541     if(BjoX>=-1)//x>0.1
542     {
543     pid_highx_parton_01_castor->Fill(i->flavour());
544     }
545     if(BjoX>=-0.698)//X>0.2
546     {
547     pid_highx_parton_02_castor->Fill(i->flavour());
548     }
549     if(BjoX>=-0.397)//x>0.4
550     {
551     pid_highx_parton_04_castor->Fill(i->flavour());
552     }
553     if(BjoX>=-0.222)//x>0.6
554     {
555     pid_highx_parton_06_castor->Fill(i->flavour());
556     }
557    
558     distX_Matching_parton_CASTOR->Fill( BjoX,fGenJet.et());
559     pid_Allx_parton_CASTOR->Fill(i->flavour());
560    
561    
562    
563    
564     }
565     /********************************************/
566     // Matching GenJet&Parton&RecJet for CASTOR //
567     /*******************************************/
568     float distance_2 = delR(GenJETEta ,CaloJetEta,GenJETPhi,CaloJetPhi);//distance for GenJet&RecJet
569     if(distance_1<threshold_1 && distance_2<threshold && cor->et()>ETCut)
570     {
571    
572     BjoX_Rec=log10((cor->pt()/sqrts)*exp((-1)*(cor->eta())));
573    
574     distX_Matching_parton_CASTOR_RecJet->Fill(BjoX_Rec,fGenJet.et());
575    
576     }
577    
578    
579     }//CASTOR eta
580    
581     }//correct RecJet
582     }//Parton
583    
584     }//GenJEt
585    
586     }
587    
588    
589     Min_DeltaR->Fill(Min_deltaR);
590    
591     /*************************************************/
592     // REC0 JET //
593     /*************************************************/
594    
595     vector<CaloJet> tmpCaloJet;
596     tmpCaloJet.clear();
597     if(CaloJetsHandle->size()){
598     for(CaloJetCollection::const_iterator it=CaloJetsHandle->begin();it!=CaloJetsHandle->end();it++)
599     tmpCaloJet.push_back(*it);
600     std::stable_sort(tmpCaloJet.begin(),tmpCaloJet.end(),CaloJetSort());
601     }
602     //cout << "There are " << tmpCaloJet.size() << " CaloJets" <<endl;
603    
604     // CaloJet TotRecJet;
605     CaloJet fRJ;
606     if(tmpCaloJet.size()>0){
607     for(int r = 0; r <tmpCaloJet.size();++r)
608     {
609     fRJ = tmpCaloJet[r];
610     if(fRJ.et()>ETCut)
611     {
612     ptCalo->Fill(fRJ.pt());
613     yCalo->Fill(fRJ.eta());
614     phiCalo->Fill(fRJ.phi());
615     ECalo->Fill(fRJ.et());
616     }
617    
618     if (abs(fRJ.eta())>=3. && abs(fRJ.eta())<=5.)
619     {
620     if(fRJ.et()>ETCut)
621     {
622    
623     ECalo_hf->Fill(fRJ.et());
624     }
625     }
626    
627     }
628     }
629    
630     /*************************************************/
631     // CORRECTED REC0 JET //
632     /*************************************************/
633     Handle<CaloJetCollection> CorJets;
634     try{
635     e.getByLabel(CorJetAlgorithm, CorJets );
636     }catch(cms::Exception& ex){
637     LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
638     }
639    
640     for(CaloJetCollection::const_iterator corjets=CorJets->begin();corjets!=CorJets->end();corjets++)
641     {
642    
643     if(corjets->et()>ETCut)
644     {
645     corr_ptCalo->Fill(corjets->pt());
646     corr_ECalo->Fill(corjets->et());
647     }
648    
649     if (abs(corjets->eta())>=3. && abs(corjets->eta())<=5.)
650     {
651     if(corjets->et()>ETCut)
652     {
653    
654     corr_ECalo_hf->Fill(corjets->et());
655     }
656     }
657     }
658    
659    
660    
661     /*******************************************/
662     // Matching dR<0.2 //
663     /******************************************/
664    
665     GenJet fGJet;
666     if(tmpGenJet.size()>0 )
667     {
668     for(int gg = 0; gg <tmpGenJet.size();++gg)
669     {
670     fGJet = tmpGenJet[gg];
671    
672     float GenJetEta=fGJet.eta();
673     float GenJetPhi=fGJet.phi();
674    
675    
676     Handle<CaloJetCollection> corjets;
677     try{
678     e.getByLabel(CorJetAlgorithm, corjets );
679     }catch(cms::Exception& ex){
680     LogError("ForwardJet") << "Error! can't get collection CaloJetsHandle";
681     }
682    
683     for(CaloJetCollection::const_iterator cor=corjets->begin();cor!=corjets->end();cor++)
684     {
685    
686     float CaloJetEta=cor->eta();
687     float CaloJetPhi=cor->phi();
688    
689     float distance = delR(GenJetEta ,CaloJetEta,GenJetPhi,CaloJetPhi);
690    
691     if (distance<=threshold && fGJet.et()>ETCut )
692     {
693    
694     RatioE_vs_Eta->Fill(abs(GenJetEta), (cor->et())/fGJet.et());
695     RatioE_vs_E->Fill(fGJet.et(), (cor->et())/fGJet.et());
696    
697     if (abs(cor->eta())>=3. && abs(cor->eta())<=5.)// FOR HF
698     {
699     double ESim= fGJet.et();
700     int trueIndex=determineIndex(limBins,ESim,nbreBins);
701     if(trueIndex!=-1)
702     {
703    
704    
705     EThist[trueIndex]->Fill((cor->et())/(fGJet.et()));
706     Phihist[trueIndex]->Fill(abs(cor->phi())/abs(fGJet.phi()));
707     Etahist[trueIndex]->Fill(abs(cor->eta())/abs(fGJet.eta()));
708     trueRecoJet->Fill(cor->et());
709     countD++;
710     }
711     RatioE_hf_vs_Eta->Fill(abs(GenJetEta), (cor->et())/fGJet.et());
712     RatioE_hf_vs_E->Fill(fGJet.et(),(cor->et()) /fGJet.et());
713    
714     }
715     }
716    
717    
718     // for not corrected calo jet
719    
720     /* for(int rr = 0; rr <tmpCaloJet.size();++rr){
721     fRJet = tmpCaloJet[rr];
722    
723     float CaloJetEta=fRJet.eta();
724     float CaloJetPhi=fRJet.phi();
725    
726     float distance = delR(GenJetEta ,CaloJetEta,GenJetPhi,CaloJetPhi);
727    
728     // if distance < threshold , fill histogramm
729     if (distance<=threshold && fGJet.et()>ETCut )
730     {
731    
732     //if(genjetused[gg] == 0)
733     //{
734    
735     RatioE_vs_Eta->Fill(abs(GenJetEta), fRJet.et()/fGJet.et());
736     RatioE_vs_E->Fill(fGJet.et(), fRJet.et()/fGJet.et());
737    
738     if (abs(fRJet.eta())>=3. && abs(fRJet.eta())<=5.)// FOR HF
739     {
740     double ESim= fGJet.et();
741     int trueIndex=determineIndex(limBins,ESim,nbreBins);
742     if(trueIndex!=-1)
743     {
744     cout<<"RjetEta=="<<fRJet.eta()<<" GenJetEt=="<<fGJet.et()<<endl;
745    
746     EThist[trueIndex]->Fill((fRJet.et())/(fGJet.et()));
747     Phihist[trueIndex]->Fill(abs(fRJet.phi())/abs(fGJet.phi()));
748     Etahist[trueIndex]->Fill(abs(fRJet.eta())/abs(fGJet.eta()));
749     trueRecoJet->Fill(fRJet.et());
750     countD++;
751     }
752     RatioE_hf_vs_Eta->Fill(abs(GenJetEta), fRJet.et()/fGJet.et());
753     RatioE_hf_vs_E->Fill(fGJet.et(), fRJet.et()/fGJet.et());
754    
755     }
756     }
757     }*/
758     }
759     }
760     }
761    
762    
763     //RecGenTree->Fill();
764     //cout << endl;
765     //cout << "===> " << e.id() << " finished" << endl;
766     eventnum++;
767     return ;
768    
769     }
770    
771     void ForwardJet::endJob(){
772 salim 1.2 long int sigma=146e7;//pb for pt_15_20 from crossection_QCD_event-list
773 salim 1.1
774     ptGen->Scale(sigma/(eventnum*ptGen->GetBinWidth(1)));
775     ptParton->Scale(sigma/(eventnum*ptParton->GetBinWidth(1)));
776     yGen->Scale(sigma/(eventnum*yGen->GetBinWidth(1)));
777     phiGen->Scale(sigma/(eventnum*phiGen->GetBinWidth(1)));
778     EGen->Scale(sigma/(eventnum* EGen->GetBinWidth(1)));
779     ptCalo->Scale(sigma/(eventnum*ptCalo->GetBinWidth(1)));
780     corr_ptCalo->Scale(sigma/(eventnum*corr_ptCalo->GetBinWidth(1)));
781     yCalo->Scale(sigma/(eventnum*yCalo->GetBinWidth(1)));
782     phiCalo->Scale(sigma/(eventnum*phiCalo->GetBinWidth(1)));
783     ECalo->Scale(sigma/(eventnum*ECalo->GetBinWidth(1)));
784     corr_ECalo->Scale(sigma/(eventnum*corr_ECalo->GetBinWidth(1)));
785     EGen_hf->Scale(sigma/(eventnum*EGen_hf->GetBinWidth(1)));
786     ECalo_hf->Scale(sigma/(eventnum*ECalo_hf->GetBinWidth(1)));
787     corr_ECalo_hf->Scale(sigma/(eventnum*corr_ECalo_hf->GetBinWidth(1)));
788     fOutputFile->Write() ;
789     fOutputFile->Close() ;
790     return ;
791     }
792     float ForwardJet::delR(const float& eta1,const float& eta2,const float& phi1,const float& phi2)
793     {
794    
795     float Distance =sqrt((2*atan(tan((phi1-phi2)/2)))*(2*atan(tan((phi1-phi2)/2)))+(eta1-eta2)*(eta1-eta2));
796     return Distance;
797    
798     }
799     int ForwardJet::determineIndex( double* limBins, double data, int nbreBins){
800    
801     int result =-1;
802     if(data<limBins[0]) return result;
803    
804     for(int i = 0; i<nbreBins;i++){
805     if(limBins[i]<= data && data<limBins[i+1]){
806     result=i;
807     break;}
808     result=i;
809     }
810    
811     return result;
812    
813     }
814    
815     DEFINE_FWK_MODULE(ForwardJet);