ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/auterman/Demo/PATJetIDAnalyzer/src/PATJetIDAnalyzer.cc
Revision: 1.1.1.2 (vendor branch)
Committed: Mon Apr 14 12:45:47 2008 UTC (17 years ago) by auterman
Content type: text/plain
Branch: tex, JetID, Demo
CVS Tags: start
Changes since 1.1.1.1: +360 -30 lines
Log Message: 
A JetID skeleton

File Contents

# User Rev Content
1 auterman 1.1 // -*- C++ -*-
2     //
3     // Package: PATJetIDAnalyzer
4     // Class: PATJetIDAnalyzer
5     //
6     /**\class PATJetIDAnalyzer PATJetIDAnalyzer.cc Demo/PATJetIDAnalyzer/src/PATJetIDAnalyzer.cc
7    
8     Description: <one line class summary>
9    
10     Implementation:
11     <Notes on implementation>
12     */
13     //
14     // Original Author: Christian Autermann
15     // Created: Mon Feb 25 11:33:02 CET 2008
16 auterman 1.1.1.2 // $Id: PATJetIDAnalyzer.cc,v 1.3 2008/04/14 08:46:50 auterman Exp $
17 auterman 1.1 //
18     //
19    
20     //system includes
21     #include <fstream>
22     #include <iostream>
23     //framework includes
24     #include "FWCore/Utilities/interface/EDMException.h"
25     #include "FWCore/ServiceRegistry/interface/Service.h"
26     #include "PhysicsTools/UtilAlgos/interface/TFileService.h"
27     //root includes
28     #include "TH1F.h"
29     #include "TH2F.h"
30     #include "TString.h"
31 auterman 1.1.1.2 #include "TRandom.h"
32 auterman 1.1 //User include files
33     #include "DataFormats/Math/interface/deltaR.h"
34     //AOD
35     #include "DataFormats/JetReco/interface/GenJet.h"
36     #include "DataFormats/JetReco/interface/GenJetfwd.h"
37     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
38     #include "DataFormats/METReco/interface/CaloMETCollection.h"
39     #include "DataFormats/METReco/interface/GenMETCollection.h"
40     #include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
41     //PAT
42     #include "DataFormats/PatCandidates/interface/Jet.h"
43     #include "DataFormats/PatCandidates/interface/MET.h"
44    
45 auterman 1.1.1.2 #include "PhysicsTools/Utilities/interface/EtComparator.h"
46 auterman 1.1 #include "Demo/PATJetIDAnalyzer/interface/NameScheme.h"
47     #include "Demo/PATJetIDAnalyzer/interface/PATJetIDAnalyzer.h"
48 auterman 1.1.1.2 #include "TMatrixT.h"
49 auterman 1.1
50 auterman 1.1.1.2 //EM-cells
51     #include "DataFormats/EcalDetId/interface/EBDetId.h"
52     #include "Geometry/Records/interface/IdealGeometryRecord.h"
53     #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
54     #include "FWCore/Framework/interface/ESHandle.h"
55     #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
56     #include "FWCore/Framework/interface/EventSetup.h"
57 auterman 1.1
58    
59 auterman 1.1.1.2 const unsigned PATJetIDAnalyzer::N_Fourier_Bins_2D;
60     const unsigned PATJetIDAnalyzer::N_Fourier_Bins_1D;
61     const unsigned PATJetIDAnalyzer::Nf_Fourier_Bins_1D;
62 auterman 1.1 //
63     // constructors and destructor
64     //
65     PATJetIDAnalyzer::PATJetIDAnalyzer(const edm::ParameterSet& iConfig) :
66 auterman 1.1.1.2 _patJet ( iConfig.getParameter<edm::InputTag>( "patJet" ) ),
67     _patMet ( iConfig.getParameter<edm::InputTag>( "patMet" ) ),
68     _ebrechits( iConfig.getParameter<edm::InputTag>("EBRecHits") ),
69     _hist ( iConfig.getParameter<std::string>( "hist" ) ),
70     _jetminpt ( iConfig.getParameter<double>( "MinJetPt" ) ),
71     _jetmaxeta(iConfig.getParameter<double>( "MaxJetEta" ) ),
72     _simulate_noise(iConfig.getParameter<bool>( "SimulateNoise" ) ),
73     _NoiseMean(iConfig.getParameter<double>( "NoiseMean" ) ),
74     _NoiseSigma(iConfig.getParameter<double>( "NoiseSigma" ) ),
75     _NoiseThreshold(iConfig.getParameter<double>( "NoiseThreshold" ) ),
76     _DoNormalization(iConfig.getParameter<bool>( "Normalize" ) ),
77     _uniqueplotid(0)
78 auterman 1.1 {
79     //now do what ever initialization is needed
80 auterman 1.1.1.2 random = new TRandom(123455);
81 auterman 1.1 }
82    
83    
84     PATJetIDAnalyzer::~PATJetIDAnalyzer()
85     {
86    
87     // do anything here that needs to be done at desctruction time
88     // (e.g. close files, deallocate resources etc.)
89    
90     }
91    
92    
93     //
94     // member functions
95     //
96    
97 auterman 1.1.1.2 double PATJetIDAnalyzer::norm(int n)
98     {
99     return ( n==0 ? 1 : sqrt(2.) );
100     }
101    
102     TH2F * PATJetIDAnalyzer::FDCT(TH2F & k)
103     {
104     TH2F* result= new TH2F(k);
105     unsigned int Nx=k.GetNbinsX();
106     unsigned int Ny=k.GetNbinsY();
107    
108     for (unsigned int u=0; u<=Nx; ++u){
109     for (unsigned int v=0; v<=Ny; ++v){
110     result->SetBinContent(u,v,0.0);
111     for (unsigned int i=0; i<Nx; ++i){
112     for (unsigned int j=0; j<Ny; ++j){
113     result->SetBinContent(u,v, result->GetBinContent(u,v)+
114     norm(u)*norm(v)/(sqrt((double)Nx)*sqrt((double)Ny))*
115     (k.GetBinContent(i,j)-0.)*
116     cos( ( (double)i+0.5)*3.14159*(double)u/( (double)Nx) )*
117     cos( ( (double)j+0.5)*3.14159*(double)v/( (double)Ny) )
118     );
119     //cout << "x="<<u<<", y="<<v<<", f="<<result->GetBinContent(u,v)<<endl;
120     }
121     }
122     }
123     }
124     return result;
125     }
126    
127     void PATJetIDAnalyzer::FakeNoise(TH2F & hk, double * k)
128     {
129     unsigned int Nx=hk.GetNbinsX();
130     unsigned int Ny=hk.GetNbinsY();
131     unsigned int r;
132     double noise, absnoise=0., abs=0.;
133     double norm = hk.Integral();
134     for (unsigned i=0; i<N_Fourier_Bins_1D; ++i)
135     abs+=k[i];
136     for (unsigned int u=1; u<=Nx; ++u){
137     for (unsigned int v=1; v<=Ny; ++v){
138     noise = random->Gaus(_NoiseMean,_NoiseSigma);
139     if (noise+hk.GetBinContent(u,v)>_NoiseThreshold) {
140     hk.SetBinContent(u,v, noise+hk.GetBinContent(u,v) );
141     double phi = ((double)(u-1)-(double)Nx)/2/(double)Nx;
142     double eta = ((double)(v-1)-(double)Ny)/2/(double)Ny;
143     r = (int)(sqrt( (phi*phi+eta*eta)/2. )*(double)(N_Fourier_Bins_1D-1)/0.5);
144     //cout <<"u="<<u<<", v="<<v<<", phi="<<phi<<", eta="<<eta<<", r="<<r<<", E(r)="<<k[r]
145     // <<", E(u,v)="<<hk.GetBinContent(u,v)<<", noise="<<noise<<endl;
146     if (r>0 && r<N_Fourier_Bins_1D){ k[r] += noise; absnoise+=noise; }
147     }
148     }
149     }
150    
151     //normalize such, that the energy of the jet ("brightness") is unchanged by noise
152     if (_DoNormalization){
153     hk.Scale( norm/hk.Integral() );
154     for (unsigned i=0; i<N_Fourier_Bins_1D; ++i)
155     k[i]=k[i]*abs/(abs+absnoise);
156     }
157     }
158    
159     void PATJetIDAnalyzer::FakeNoise2x2(TH2F & hk, double * k)
160     {
161     unsigned int Nx=hk.GetNbinsX();
162     unsigned int Ny=hk.GetNbinsY();
163     unsigned int r;
164     double noise, absnoise=0., abs=0.;
165     double norm = hk.Integral();
166     for (unsigned i=0; i<N_Fourier_Bins_1D; ++i)
167     abs+=k[i];
168     unsigned int u=(int)(random->Rndm()*(double)hk.GetNbinsX())+1;
169     unsigned int v=(int)(random->Rndm()*(double)hk.GetNbinsY())+1;
170     noise = 15.;
171     cout << u << ", " << v << ": old="<<hk.GetBinContent(u,v);
172     hk.SetBinContent(u,v, noise+hk.GetBinContent(u,v) );
173     cout << "; new="<<hk.GetBinContent(u,v)<<endl;;
174     double phi = ((double)(u-1)-(double)Nx)/2/(double)Nx;
175     double eta = ((double)(v-1)-(double)Ny)/2/(double)Ny;
176     r = (int)(sqrt( (phi*phi+eta*eta)/2. )*(double)(N_Fourier_Bins_1D-1)/0.5);
177     //cout <<"u="<<u<<", v="<<v<<", phi="<<phi<<", eta="<<eta<<", r="<<r<<", E(r)="<<k[r]
178     // <<", E(u,v)="<<hk.GetBinContent(u,v)<<", noise="<<noise<<endl;
179     if (r>0 && r<N_Fourier_Bins_1D){ k[r] += noise; absnoise+=noise; }
180    
181     //normalize such, that the energy of the jet ("brightness") is unchanged by noise
182     if (_DoNormalization){
183     hk.Scale( norm/hk.Integral() );
184     for (unsigned i=0; i<N_Fourier_Bins_1D; ++i)
185     k[i]=k[i]*abs/(abs+absnoise);
186     }
187     }
188    
189     double PATJetIDAnalyzer::dphi(double phi1, double phi2)
190     {
191     double result = phi1 - phi2;
192     if (result> 3.14159) result-=3.14159;
193     if (result<-3.14159) result+=3.14159;
194     return result;
195     }
196    
197    
198     void PATJetIDAnalyzer::FourierTransformation( const unsigned int i,
199     const pat::Jet& jet,
200     const EBRecHitCollection& EBRecHit)
201     {
202     ///Studies
203     double k[N_Fourier_Bins_1D], f[Nf_Fourier_Bins_1D];
204     for (unsigned b=0; b<N_Fourier_Bins_1D; ++b){
205     k[b]=0.;
206     //k[b] = cos(2.*3.1415*(double)b/(double)N_Fourier_Bins) +
207     // cos(0.2*2.*3.1415*(double)b/(double)N_Fourier_Bins) +
208     // cos(10*2.*3.1415*(double)b/(double)N_Fourier_Bins);
209    
210     }
211     for (unsigned b=0; b<Nf_Fourier_Bins_1D; ++b){
212     f[b]=0.;
213     }
214     char * name = new char[64];
215     sprintf(name,"h_%d",++_uniqueplotid);
216     TH2F * k2d = new TH2F(name,"",16,-8,7,16,-8,7);
217     TH2F * f2d;
218    
219     int jtow=0; //, icell=0;
220     std::vector <CaloTowerRef> jetTowers = jet.getConstituents();
221    
222     for(std::vector<CaloTowerRef>::const_iterator tow = jetTowers.begin();
223     tow != jetTowers.end(); ++tow, ++jtow){
224    
225     //2D-energy picture of this jet
226     k2d->Fill( dphi((*tow)->phi(),jet.phi())/0.0873,
227     ((*tow)->eta()-jet.eta())/0.087,
228     (*tow)->energy() );
229    
230     //1D-energy picture
231     double dR = deltaR(jet,*(*tow));
232     if (fabs(dR)>0.5) dR=0.5;
233     //_ft_k[i]->Fill( dR, (*tow)->energy() );
234     k[ (int)(((double)N_Fourier_Bins_1D-1.)*dR/0.5) ] += (*tow)->energy();
235    
236     //no EM cells in my current test data....
237     /*
238     double eem=0.;
239     for (size_t it=0; it<(*tow)->constituentsSize(); ++it) {
240     const DetId detid = (*tow)->constituent(it);
241     EcalRecHitCollection::const_iterator myRecHit = EBRecHit.find(detid);
242     if(myRecHit != EBRecHit.end()) {
243     eem += myRecHit->energy();
244     EBDetId det = myRecHit->id();
245    
246     const CaloCellGeometry* cell=EBgeom->getGeometry( myRecHit->id() );
247     /
248     etowet [icell] = myRecHit->energy()*sin( cell->getPosition().theta());
249     etoweta[icell] = cell->getPosition().eta();
250     etowphi[icell] = cell->getPosition().phi();
251     etowe [icell] = myRecHit->energy();
252     etowid_phi[icell] = det.iphi();
253     etowid_eta[icell] = det.ieta();
254     etowid [icell] = myRecHit->id().rawId();
255     etownum[icell] = icell;
256     //
257     _ft_energy[i]->Fill( deltaPhi( (double)cell->getPosition().phi(),jet.phi()), cell->getPosition().eta()-jet.eta(), myRecHit->energy() );
258     ++icell;
259     }
260     }
261     */
262     }
263    
264    
265     //Add fake-noise to the jet:
266     if(_simulate_noise)
267     FakeNoise( *k2d, k);
268     //FakeNoise2x2( *k2d, k);
269    
270    
271     //1-D fourier transformation
272     for (unsigned b=0; b<N_Fourier_Bins_1D; ++b) {
273     _ft_k[i]->SetBinContent(b, k[b] + _ft_k[i]->GetBinContent(b) );
274     //before fourier-trafo substract average jet spectrum:
275     k[b] -= 167.53*exp( -11.941*0.5*(double)b/(double)N_Fourier_Bins_1D );
276     _ft_ksubavg[i]->SetBinContent(b, k[b] + _ft_ksubavg[i]->GetBinContent(b) );
277     }
278     for (unsigned b=0; b<=Nf_Fourier_Bins_1D; ++b) {
279     double fourier=0.;
280     for (unsigned j=0; j<N_Fourier_Bins_1D; ++j) {
281     fourier += 2./0.5*k[j]*
282     cos( ((double)j)*3.14159*(double)b/
283     (double)N_Fourier_Bins_1D );
284     }
285     //_ft_f[i]->SetBinContent(b, fourier);
286     _ft_f[i]->SetBinContent(b+1, fourier + _ft_f[i]->GetBinContent(b+1) );
287     }
288    
289    
290     //2D-fourier transformation & calculation of observables
291     TMatrixT<double> matrix16(16,16);
292     TMatrixT<double> matrix8(8,8);
293     for (int binx=0; binx<=k2d->GetNbinsX();++binx) {
294     for (int biny=0; biny<=k2d->GetNbinsY();++biny) {
295     _ft_energy[i]->SetBinContent(binx,biny,_ft_energy[i]->GetBinContent(binx,biny)+
296     k2d->GetBinContent(binx,biny) );
297     }}
298    
299     std::vector<double> db;
300     double sum=0.;
301     f2d = FDCT( *k2d);
302     for (int binx=0; binx<=f2d->GetNbinsX();++binx) {
303     for (int biny=0; biny<=f2d->GetNbinsY();++biny) {
304     double bincont = f2d->GetBinContent(binx,biny);
305     _ft_frequency[i]->SetBinContent(binx,biny, _ft_frequency[i]->GetBinContent(binx,biny)+
306     bincont );
307     if (bincont!=0.) db.push_back( fabs(bincont) );
308     sum += fabs(bincont);
309     if (binx>0 && binx<16 &&biny>0 && biny<16)
310     matrix16[binx][biny] = f2d->GetBinContent(binx,biny);
311     if (binx>3 && binx<12 &&biny>3 && biny<12) {
312     matrix8[binx-4][biny-4] = f2d->GetBinContent(binx,biny);
313     //cout <<"(" <<binx << "::"<<biny<< ") = " <<k2d->GetBinContent(bin)<<endl;
314     }
315     }}
316     std::sort(db.begin(), db.end());
317     double ig=0., F10=0., Fs10=0.;
318     int n10=-1, n30=-1, n60=-1, n90=-1, n95=-1, n99=-1;
319     for (unsigned j=0; j<db.size(); ++j){
320     if ( ig>0.1*sum && n10==-1) n10 = j;
321     if ( ig>0.3*sum && n30==-1) n30 = j;
322     if ( ig>0.6*sum && n60==-1) n60 = j;
323     if ( ig>0.9*sum && n90==-1) n90 = j;
324     if ( ig>0.95*sum && n95==-1) n95 = j;
325     if ( ig>0.99*sum && n99==-1) n99 = j;
326     ig += fabs(db[j]);
327     if (j<10) Fs10 += db[j];
328     if (j>=db.size()-10 && db.size()>10) F10 += db[j];
329     }
330     double Ilow = f2d->Integral( 1, (int)(f2d->GetNbinsX()/3.),0,(int)(f2d->GetNbinsY()/3.) );
331     double Ihi = f2d->Integral( (int)(f2d->GetNbinsX()*2./3.),f2d->GetNbinsX(), (int)(f2d->GetNbinsY()*2./3.), f2d->GetNbinsY()-1 );
332    
333     //fill observables from the 2D-fourier-transformation:
334     _fto_n99[i]->Fill( n99 );
335     _fto_n95[i]->Fill( n95 );
336     _fto_n90[i]->Fill( n90 );
337     _fto_n60[i]->Fill( n60 );
338     _fto_n30[i]->Fill( n30 );
339     _fto_n10[i]->Fill( n10 );
340     _fto_F10[i]->Fill( F10 );
341     _fto_Fs10[i]->Fill( Fs10 );
342     if (Ihi!=0.) {
343     _fto_LowFovHiF[i]->Fill( Ilow/Ihi );
344     _fto_LowFvsHiF[i]->Fill( (Ihi-Ilow)/(Ihi+Ilow) );
345     }
346     //_fto_det16[i]->Fill( matrix16.Determinant() );
347     //_fto_det8[i]->Fill( matrix8.Determinant() );
348    
349     //clean-up
350     delete name;
351     delete f2d;
352     delete k2d;
353     }
354    
355 auterman 1.1 // ------------ method called to for each event ------------
356     void
357     PATJetIDAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
358     {
359     using namespace edm;
360 auterman 1.1.1.2 using namespace reco;
361 auterman 1.1 using namespace std;
362 auterman 1.1.1.2 using namespace pat;
363    
364     // edm::ESHandle<CaloGeometry> pG;
365     // iSetup.get<IdealGeometryRecord>().get(pG);
366     // const CaloGeometry cG = *pG;
367     // EBgeom=cG.getSubdetectorGeometry(DetId::Ecal,1);
368    
369     //EB RecHits
370     const EBRecHitCollection *EBRecHit = 0;
371     edm::Handle<EBRecHitCollection> EcalRecHitEB;
372     // iEvent.getByLabel( _ebrechits, EcalRecHitEB);
373     // if( EcalRecHitEB.isValid() ){
374     // EBRecHit = EcalRecHitEB.product();
375     // }
376    
377     //PAT Jets
378     typedef vector<pat::Jet>::const_iterator PatJetCIter;
379     edm::Handle<vector<pat::Jet> > PatJets;
380     iEvent.getByLabel( _patJet, PatJets );
381     GreaterByEt<pat::Jet> eTComparator_;
382     vector<pat::Jet> patjets = *PatJets;
383     std::sort(patjets.begin(), patjets.end(), eTComparator_);
384     PatJetCIter PatJetBegin = patjets.begin();
385     PatJetCIter PatJetEnd = patjets.end();
386 auterman 1.1
387 auterman 1.1.1.2 //Jets
388 auterman 1.1 unsigned multiplicity = 0;
389     double met=0.0, ht=0.0;
390 auterman 1.1.1.2 for (PatJetCIter it=PatJetBegin; it!=PatJetEnd; ++it) {
391 auterman 1.1 if (it->pt() > _jetminpt && fabs( it->eta() )<_jetmaxeta){
392     if (multiplicity<_njets ) {
393     _pt_jet[ multiplicity]->Fill( it->pt() );
394     _eta_jet[ multiplicity]->Fill( it->eta() );
395     _phi_jet[ multiplicity]->Fill( it->phi() );
396     _emfrac_jet[ multiplicity]->Fill( it->emEnergyFraction() );
397     _hadfrac_jet[multiplicity]->Fill( it->energyFractionHadronic() );
398 auterman 1.1.1.2
399 auterman 1.1 _n60_jet[ multiplicity]->Fill( it->n60() );
400     _n90_jet[ multiplicity]->Fill( it->n90() );
401     _area_jet[ multiplicity]->Fill( it->towersArea() );
402 auterman 1.1.1.2
403     FourierTransformation( multiplicity, *it, *EBRecHit );
404 auterman 1.1 }
405     ++multiplicity;
406     ht += it->pt();
407 auterman 1.1.1.2 met += it->pt();//@@ just to kill warnings....
408 auterman 1.1 }
409     }
410     _jetmult->Fill( multiplicity );
411     _ht->Fill( ht );
412 auterman 1.1.1.2 if (multiplicity>1)
413     _dijet->Fill( sqrt( pow(patjets[0].energy()+patjets[1].energy(),2) -
414     pow(patjets[0].px()+patjets[1].px(),2) -
415     pow(patjets[0].py()+patjets[1].py(),2) -
416     pow(patjets[0].pz()+patjets[1].pz(),2)
417     ) );
418 auterman 1.1
419 auterman 1.1.1.2
420    
421     /*
422     //PAT MET
423     typedef vector<PatMet>::const_iterator PatJetCIter;
424     edm::Handle<vector<PatJet> > PatJets;
425     iEvent.getByLabel( _patJet, PatJets );
426     std::sort( PatJets->begin(), PatJets->end(), PtGreater());
427     PatJetCIter PatJetBegin = PatJets->begin();
428     PatJetCIter PatJetEnd = PatJets->end();
429 auterman 1.1 //MET
430     edm::Handle<reco::CaloMETCollection> recMet;
431     iEvent.getByLabel(_recMet,recMet);
432    
433     for ( reco::CaloMETCollection::const_iterator it=recMet->begin();
434     it!=recMet->end(); ++it) {
435     met = it->pt();
436     _met->Fill( met );
437     _metx->Fill(met*sin(it->phi()));
438     _mety->Fill(met*cos(it->phi()));
439     break;
440     }
441     _metmult->Fill( recMet->size() );
442     _htmet->Fill( ht+met );
443     _metsig->Fill( met / (0.97*sqrt(ht)) );
444     if (multiplicity>1 ) {
445     _nv->Fill( (met*met-(CaloJets[0].pt()-CaloJets[1].pt())*(CaloJets[0].pt()-CaloJets[1].pt())) / (CaloJets[1].pt()*CaloJets[1].pt()) );
446     _nv2->Fill( (met*met-(CaloJets[0].pt()-CaloJets[1].pt())*(CaloJets[0].pt()-CaloJets[1].pt())) / (CaloJets[0].pt()*CaloJets[0].pt()) );
447     _dijet->Fill(sqrt(pow(CaloJets[0].energy()+CaloJets[1].energy(),2)
448     -(pow(CaloJets[0].px()+CaloJets[1].px(),2)
449     +pow(CaloJets[0].py()+CaloJets[1].py(),2)
450     +pow(CaloJets[0].pz()+CaloJets[1].pz(),2) ) ) );
451     }
452 auterman 1.1.1.2 */
453     /*
454 auterman 1.1 //GenJets
455     multiplicity = 0;
456     met=0.0;
457     ht=0.0;
458     for (reco::GenJetCollection::const_iterator it=GenJets.begin();
459     it!=GenJets.end(); ++it) {
460     if (it->pt() > _jetminpt && fabs( it->eta() )<_jetmaxeta) {
461     if (multiplicity<_ngenjets ) {
462     _pt_genjet[ multiplicity]->Fill( it->pt() );
463     _eta_genjet[ multiplicity]->Fill( it->eta() );
464     _phi_genjet[ multiplicity]->Fill( it->phi() );
465     _emfrac_genjet[ multiplicity]->Fill( it->emEnergy() / it->energy() );
466     _hadfrac_genjet[ multiplicity]->Fill( it->hadEnergy() / it->energy() );
467     _invisible_genjet[multiplicity]->Fill( it->invisibleEnergy() / it->energy() );
468     _aux_genjet[ multiplicity]->Fill( it->auxiliaryEnergy() / it->energy() );
469     }
470     ++multiplicity;
471     ht += it->pt();
472     }
473     }
474     _genjetmult->Fill( multiplicity );
475    
476     //GenMET
477     edm::Handle<reco::GenMETCollection> genMet;
478     iEvent.getByLabel(_genMet,genMet);
479    
480     for ( reco::GenMETCollection::const_iterator it=genMet->begin();
481     it!=genMet->end(); ++it) {
482     met = it->pt();
483     _genmet->Fill( met );
484     _genmetx->Fill(met*sin(it->phi()));
485     _genmety->Fill(met*cos(it->phi()));
486     break;
487     }
488     _genmetmult->Fill( recMet->size() );
489     _genht->Fill( ht );
490     _genhtmet->Fill( ht+met );
491     if (multiplicity>1 ) {
492     _gendijet->Fill(sqrt(pow(GenJets[0].energy()+GenJets[1].energy(),2)
493     -(pow(GenJets[0].px()+GenJets[1].px(),2)
494     +pow(GenJets[0].py()+GenJets[1].py(),2)
495     +pow(GenJets[0].pz()+GenJets[1].pz(),2) ) ) );
496     }
497 auterman 1.1.1.2 */
498    
499    
500 auterman 1.1
501     }
502    
503    
504     // ------------ method called once each job just before starting event loop ------------
505     void
506     PATJetIDAnalyzer::beginJob(const edm::EventSetup&)
507     {
508     if( _hist.empty() )
509     return;
510    
511     edm::Service<TFileService> fs;
512     if( !fs ){
513     throw edm::Exception( edm::errors::Configuration,
514     "TFile Service is not registered in cfg file" );
515     }
516    
517 auterman 1.1.1.2 NameScheme hsusy("pat");
518 auterman 1.1 ofstream hist(_hist.c_str(), std::ios::out);
519     //CaloJets
520     for (unsigned i=0; i<_njets; ++i) {
521     _pt_jet[i] = fs->make<TH1F>(hsusy.name(hist, "ptj", i ), hsusy.name("PtJet", i), 100, 0. , 1000.);
522     _eta_jet[i] = fs->make<TH1F>(hsusy.name(hist, "etaj", i ), hsusy.name("EtaJet", i), 100, -5. , 5.);
523     _phi_jet[i] = fs->make<TH1F>(hsusy.name(hist, "phij", i ), hsusy.name("PhiJet", i), 70, -3.5, 3.5);
524     _emfrac_jet[i] = fs->make<TH1F>(hsusy.name(hist, "emfj", i ), hsusy.name("EmfJet", i), 50, 0. , 2.);
525     _hadfrac_jet[i] = fs->make<TH1F>(hsusy.name(hist, "hadfj",i ), hsusy.name("HadfJet",i), 50, 0. , 2.);
526     _n60_jet[i] = fs->make<TH1F>(hsusy.name(hist, "n60j", i ), hsusy.name("n60Jet", i), 50, 0. , 50.);
527     _n90_jet[i] = fs->make<TH1F>(hsusy.name(hist, "n90j", i ), hsusy.name("n90Jet", i), 50, 0. , 50.);
528     _area_jet[i] = fs->make<TH1F>(hsusy.name(hist, "areaj",i ), hsusy.name("AreaJet",i), 50, 0. , 5.);
529 auterman 1.1.1.2
530     _ft_energy[i] = fs->make<TH2F>(hsusy.name(hist, "ft_E_j",i ), hsusy.name("FtEJet",i), 16,-8,8,16,-8,8);
531     _ft_frequency[i]= fs->make<TH2F>(hsusy.name(hist, "ft_F_j",i ), hsusy.name("FtFJet",i), 16,-8,8,16,-8,8);
532     //_ft_frequency[i]= fs->make<TH1F>(hsusy.name(hist, "ft_F_j",i ), hsusy.name("FtFJet",i), 100,-50,50);
533     _ft_k[i] = fs->make<TH1F>(hsusy.name(hist, "ft_k_j",i ), hsusy.name("FtKJet",i), N_Fourier_Bins_1D,0.,0.5);
534     _ft_f[i] = fs->make<TH1F>(hsusy.name(hist, "ft_f_j",i ), hsusy.name("FtFJet",i), Nf_Fourier_Bins_1D+2,0.,N_Fourier_Bins_1D+1);
535     _ft_ksubavg[i] = fs->make<TH1F>(hsusy.name(hist, "ft_ksubavg_j",i ), hsusy.name("FtKsubavgJet",i), N_Fourier_Bins_1D,0.,0.5);
536    
537     _noisecontrib[i]= fs->make<TH1F>(hsusy.name(hist, "noisecontrib_j",i ), hsusy.name("noisecontrib",i), 100,0.,20.);
538     _fto_n99[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n99_j",i ), hsusy.name("fto_n99",i), 100,-1,500);
539     _fto_n95[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n95_j",i ), hsusy.name("fto_n95",i), 100,-1,500);
540     _fto_n90[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n90_j",i ), hsusy.name("fto_n90",i), 100,-1,500);
541     _fto_n60[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n60_j",i ), hsusy.name("fto_n60",i), 100,-1,500);
542     _fto_n30[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n30_j",i ), hsusy.name("fto_n30",i), 50,-1,500);
543     _fto_n10[i] = fs->make<TH1I>(hsusy.name(hist, "fto_n10_j",i ), hsusy.name("fto_n10",i), 50,-1,500);
544     _fto_F10[i] = fs->make<TH1F>(hsusy.name(hist, "fto_F10_j",i ), hsusy.name("fto_F10",i), 100,0.,2000.);
545     _fto_Fs10[i] = fs->make<TH1F>(hsusy.name(hist, "fto_Fs10_j",i),hsusy.name("fto_Fs10",i), 100,0.,100.);
546     _fto_LowFvsHiF[i]=fs->make<TH1F>(hsusy.name(hist, "fto_LowFvsHiF_j",i ), hsusy.name("LowFvsHiF",i), 100,-5.,5.);
547     _fto_LowFovHiF[i]=fs->make<TH1F>(hsusy.name(hist, "fto_LowFovHiF_j",i ), hsusy.name("LowFovHiF",i), 100,-5.,5.);
548     _fto_det16[i] =fs->make<TH1F>(hsusy.name(hist, "fto_Det16_j",i ), hsusy.name("Det16",i), 100,-1.0E08,1.0E08);
549     _fto_det8[i] =fs->make<TH1F>(hsusy.name(hist, "fto_Det8_j",i ), hsusy.name("Det8",i), 100,-1.0E08,1.0E08);
550    
551 auterman 1.1 }
552     _jetmult = fs->make<TH1F>(hsusy.name(hist, "multj" ), hsusy.name("JetMultiplicity"), 21, -0.5, 20.5);
553     _ht = fs->make<TH1F>(hsusy.name(hist, "ht" ), hsusy.name("Ht"), 100, 0.0, 2000.);
554     _htmet = fs->make<TH1F>(hsusy.name(hist, "htmet" ), hsusy.name("HtMet"), 100, 0.0, 2000.);
555     _dijet = fs->make<TH1F>(hsusy.name(hist, "dijet" ), hsusy.name("DiJet"), 100, 0.0, 2000.);
556     //GenJets
557     for (unsigned i=0; i<_ngenjets; ++i) {
558     _pt_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "ptgj", i ), hsusy.name("PtGenJet", i), 100, 0. , 1000.);
559     _eta_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "etagj", i ), hsusy.name("EtaGenJet", i), 100, -5. , 5.);
560     _phi_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "phigj", i ), hsusy.name("PhiGenJet", i), 70, -3.5, 3.5);
561     _emfrac_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "emfgj", i ), hsusy.name("EmfGenJet", i), 50, 0. , 5.);
562     _hadfrac_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "hadfgj",i ), hsusy.name("HadfGenJet",i), 50, 0. , 5.);
563     _invisible_genjet[i]=fs->make<TH1F>(hsusy.name(hist, "invisgj",i), hsusy.name("InvisiblefGenJet", i), 50, 0.,5.);
564     _aux_genjet[i] = fs->make<TH1F>(hsusy.name(hist, "auxgj", i ), hsusy.name("AuxfGenJet", i), 50, 0. , 5.);
565     }
566     _genjetmult = fs->make<TH1F>(hsusy.name(hist, "multgj" ), hsusy.name("GenJetMultiplicity"), 21, -0.5, 20.5);
567     _genht = fs->make<TH1F>(hsusy.name(hist, "htgj" ), hsusy.name("GenHt"), 100, 0.0, 2000.);
568     _genhtmet = fs->make<TH1F>(hsusy.name(hist, "htmetgj"), hsusy.name("GenHtMet"), 100, 0.0, 2000.);
569     _gendijet = fs->make<TH1F>(hsusy.name(hist, "dijetgj"), hsusy.name("GenDiJet"), 100, 0.0, 2000.);
570     //MET
571     _met = fs->make<TH1F>(hsusy.name(hist, "met" ), hsusy.name("Met"), 100, 0.0, 1000.);
572     _metmult = fs->make<TH1F>(hsusy.name(hist, "metmult"),hsusy.name("MetMult"), 11, -0.5, 10.5);
573     _metx = fs->make<TH1F>(hsusy.name(hist, "metx" ), hsusy.name("MetX"), 100, -1000.0, 1000.);
574     _mety = fs->make<TH1F>(hsusy.name(hist, "mety" ), hsusy.name("MetY"), 100, -1000.0, 1000.);
575     _metsig = fs->make<TH1F>(hsusy.name(hist, "metsig"), hsusy.name("MetSig"), 100, 0.0, 100.);
576     //GenMET
577     _genmet = fs->make<TH1F>(hsusy.name(hist, "metgj" ), hsusy.name("GenMet"), 100, 0.0, 1000.);
578     _genmetmult = fs->make<TH1F>(hsusy.name(hist, "metmultgj"),hsusy.name("GenMetMult"), 11, -0.5, 10.5);
579     _genmetx = fs->make<TH1F>(hsusy.name(hist, "metxgj" ), hsusy.name("GenMetX"), 100, -1000.0, 1000.);
580     _genmety = fs->make<TH1F>(hsusy.name(hist, "metygj" ), hsusy.name("GenMetY"), 100, -1000.0, 1000.);
581     //other
582     _nv = fs->make<TH1F>(hsusy.name(hist, "nv" ), hsusy.name("NV"), 100, -5.0, 5.);
583     _nv2 = fs->make<TH1F>(hsusy.name(hist, "nv2" ), hsusy.name("NV2"), 100, -5.0, 5.);
584     //matching
585     _GenOnCalo_match = fs->make<TH2F>(hsusy.name(hist,"gencalo"), hsusy.name("GenOnCalo"), 60, 0., 3., 120, -1, 5.);
586     _CaloOnGen_match = fs->make<TH2F>(hsusy.name(hist,"calogen"), hsusy.name("CaloOnGen"), 60, 0., 3., 120, -1, 5.);
587     _GenVsMatched_pt = fs->make<TH2F>(hsusy.name(hist,"genmatched_pt"), hsusy.name("GenVsMatched_Pt") , 100, 0. , 2000. , 100, 0. , 2000.);
588     _GenVsMatched_eta = fs->make<TH2F>(hsusy.name(hist,"genmatched_eta"), hsusy.name("GenVsMatched_Eta"), 100, -5. , 5. , 100, -5. , 5.);
589     _GenVsMatched_phi = fs->make<TH2F>(hsusy.name(hist,"genmatched_phi"), hsusy.name("GenVsMatched_Phi"), 70, -3.5, 3.5, 70, -3.5, 3.5);
590     _RecoEff_pt = fs->make<TH1F>(hsusy.name(hist, "recoeffpt") , hsusy.name("RecoEff_Pt") , 100, 0. , 2000. );
591     _RecoEff_eta = fs->make<TH1F>(hsusy.name(hist, "recoeffeta"), hsusy.name("RecoEff_Eta"), 100, -5. , 5. );
592     _RecoEff_phi = fs->make<TH1F>(hsusy.name(hist, "recoeffphi"), hsusy.name("RecoEff_Phi"), 70, -3.5, 3.5);
593     //helper histograms
594     _RecoTot_pt = new TH1F("recototpt" , "RecoTot_Pt" , 100, 0. , 2000. );
595     _RecoTot_eta= new TH1F("recototeta", "RecoTot_Eta", 100, -5. , 5. );
596     _RecoTot_phi= new TH1F("recototphi", "RecoTot_Phi", 70, -3.5, 3.5);
597     }
598    
599     // ------------ method called once each job just after ending the event loop ------------
600     void
601     PATJetIDAnalyzer::endJob() {
602     _RecoTot_pt->Sumw2();
603     _RecoEff_pt->Sumw2();
604     _RecoEff_pt->Divide(_RecoTot_pt);
605     _RecoTot_eta->Sumw2();
606     _RecoEff_eta->Sumw2();
607     _RecoEff_eta->Divide(_RecoTot_eta);
608     _RecoTot_phi->Sumw2();
609     _RecoEff_phi->Sumw2();
610     _RecoEff_phi->Divide(_RecoTot_phi);
611     }
612    
613     // ------------ redo matching
614     void
615     PATJetIDAnalyzer::makeMatchingMaps(edm::Handle<reco::GenJetCollection> GenJets, edm::Handle<reco::CaloJetCollection> CaloJets)
616     {
617     MatchingMapGen.clear();
618     MatchingMapJet.clear();
619    
620     // Loop over generator Jets
621     for(reco::GenJetCollection::const_iterator genjet = GenJets->begin(); genjet!=GenJets->end(); ++genjet) {
622     const reco::CaloJet* jetbest=0;
623     const reco::CaloJet* jetmatched=0;
624     double minDelR= 0.3;
625     double minDelE=-1/4.;
626     double maxDelE= 3.0;
627     double mindE=999.;
628     double mindR=999.;
629     double delR;
630     double delE;
631     double E1=(*genjet).energy();
632     bool matched = false;
633     // Loop over Calo jets
634     for(reco::CaloJetCollection::const_iterator calojet = CaloJets->begin(); calojet!=CaloJets->end(); ++calojet){
635     double E2=(*calojet).energy();
636     delR=deltaR(*calojet,*genjet);
637     delE=(E2-E1)/E1;
638     // This is the matching condition
639     if(delR<minDelR && delE>minDelE && delE<maxDelE){
640     matched=true;
641     minDelR=delR;
642     jetmatched=&(*calojet);
643     }
644     if(delR<mindR){
645     mindR=delR;
646     mindE=delE;
647     jetbest=&(*calojet);
648     }
649     }
650     if (matched) {
651     MatchingMapGen.insert(make_pair(&(*genjet),&(*jetmatched)));
652     _GenVsMatched_pt->Fill((*genjet).pt(),(*jetmatched).pt());
653     _GenVsMatched_eta->Fill((*genjet).eta(),(*jetmatched).eta());
654     _GenVsMatched_phi->Fill((*genjet).phi(),(*jetmatched).phi());
655     _RecoEff_pt->Fill((*genjet).pt());
656     _RecoEff_eta->Fill((*genjet).eta());
657     _RecoEff_phi->Fill((*genjet).phi());
658     }
659     _RecoTot_pt->Fill((*genjet).pt());
660     _RecoTot_eta->Fill((*genjet).eta());
661     _RecoTot_phi->Fill((*genjet).phi());
662     _GenOnCalo_match->Fill(mindR,mindE);
663     }
664    
665     // Loop over Calo Jets
666     for(reco::CaloJetCollection::const_iterator calojet = CaloJets->begin(); calojet!=CaloJets->end(); ++calojet) {
667     const reco::GenJet* jetbest=0;
668     const reco::GenJet* jetmatched=0;
669     double minDelR= 0.3;
670     double minDelE=-1/4.;
671     double maxDelE= 3.0;
672     double mindE=999.;
673     double mindR=999.;
674     double delR;
675     double delE;
676     double E1=(*calojet).energy();
677     bool matched = false;
678     // Loop over jets
679     for(reco::GenJetCollection::const_iterator genjet = GenJets->begin(); genjet!=GenJets->end(); ++genjet){
680     double E2=(*genjet).energy();
681     delR=deltaR(*genjet,*calojet);
682     delE=(E2-E1)/E1;
683     // This is the matching condition
684     if(delR<minDelR && delE>minDelE && delE<maxDelE){
685     matched=true;
686     minDelR=delR;
687     jetmatched=&(*genjet);
688     }
689     if(delR<mindR){
690     mindR=delR;
691     mindE=delE;
692     jetbest=&(*genjet);
693     }
694     }
695     if (matched) {
696     MatchingMapJet.insert(make_pair(&(*calojet),&(*jetmatched)));
697     }
698     _CaloOnGen_match->Fill(mindR,mindE);
699     }
700    
701     }
702    
703     //define this as a plug-in
704     DEFINE_FWK_MODULE(PATJetIDAnalyzer);