ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/AFB/src/AFB.cc
Revision: 1.6
Committed: Wed Jan 5 16:30:07 2011 UTC (14 years, 3 months ago) by efe
Content type: text/plain
Branch: MAIN
Changes since 1.5: +64 -5 lines
Log Message:
added dimuon trigger and matching

File Contents

# User Rev Content
1 efe 1.1 // -*- C++ -*-
2     //
3     // Package: AFB
4     // Class: AFB
5     //
6     /**\class AFB AFB.cc ZJet/AFB/src/AFB.cc
7    
8     Description: <one line class summary>
9    
10     Implementation:
11     <Notes on implementation>
12     */
13     //
14     // Original Author: Efe Yazgan
15     // Created: Tue Feb 3 10:08:43 CET 2009
16 efe 1.6 // $Id: AFB.cc,v 1.5 2010/11/19 16:21:10 efe Exp $
17 efe 1.1 //
18     //
19     #include <memory>
20     #include "FWCore/Framework/interface/Frameworkfwd.h"
21     #include "FWCore/Framework/interface/EDAnalyzer.h"
22     #include "FWCore/Framework/interface/Event.h"
23     #include "FWCore/Framework/interface/MakerMacros.h"
24     #include "FWCore/ParameterSet/interface/ParameterSet.h"
25     #include "FWCore/Framework/interface/Event.h"
26 efe 1.2
27     #include "FWCore/Framework/interface/EventSetupRecord.h"
28    
29 efe 1.1 #include "FWCore/Framework/interface/Run.h"
30     #include "TH1.h"
31     #include "TH1F.h"
32     #include "TH2F.h"
33     #include "TProfile.h"
34     #include <TROOT.h>
35     #include "TF1.h"
36     #include "TMath.h"
37     #include <TSystem.h>
38     #include "TFile.h"
39     #include <TCanvas.h>
40     #include <cmath>
41     #include <iostream>
42     #include <fstream>
43     #include <vector>
44     #include <functional>
45     #include <Math/VectorUtil.h>
46     #include "DataFormats/Common/interface/Handle.h"
47     #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
48     #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
49     #include "DataFormats/Candidate/interface/CompositePtrCandidate.h"
50     #include "DataFormats/Candidate/interface/Candidate.h"
51     #include "DataFormats/Candidate/interface/CandAssociation.h"
52     #include "FWCore/ServiceRegistry/interface/Service.h"
53     #include "CommonTools/UtilAlgos/interface/TFileService.h"
54     //muons
55     #include "DataFormats/MuonReco/interface/MuonFwd.h"
56     #include "DataFormats/MuonReco/interface/Muon.h"
57     #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
58     #include "DataFormats/MuonReco/interface/MuonChamberMatch.h"
59     #include "DataFormats/MuonReco/interface/MuonIsolation.h"
60     #include "DataFormats/MuonReco/interface/MuonEnergy.h"
61     #include "DataFormats/MuonReco/interface/MuonTime.h"
62     #include "DataFormats/MuonReco/interface/MuonQuality.h"
63     #include "DataFormats/MuonReco/interface/MuonSelectors.h"
64     #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
65     #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
66     #include "DataFormats/TrackReco/interface/Track.h"
67     #include "DataFormats/TrackReco/interface/TrackFwd.h"
68     #include "DataFormats/Common/interface/ValueMap.h"
69     //jets
70     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
71     #include "DataFormats/JetReco/interface/CaloJet.h"
72     #include "DataFormats/JetReco/interface/GenJetCollection.h"
73     #include "DataFormats/JetReco/interface/GenJet.h"
74     #include "JetMETCorrections/Objects/interface/JetCorrector.h"
75     //#include "RecoJets/JetAlgorithms/interface/JetIDHelper.h"
76     #include "RecoJets/JetProducers/interface/JetIDHelper.h"
77     // MET
78     #include "DataFormats/METReco/interface/CaloMET.h"
79     #include "DataFormats/METReco/interface/CaloMETFwd.h"
80     #include "DataFormats/METReco/interface/MET.h"
81     #include "DataFormats/METReco/interface/METFwd.h"
82     #include "DataFormats/METReco/interface/PFMET.h"
83     #include "DataFormats/METReco/interface/PFMETFwd.h"
84     //electron isolation
85     #include "FWCore/Utilities/interface/InputTag.h"
86     //#include "SHarper/EgIsolExample/interface/EgIsolExample.h"
87     #include "DataFormats/BeamSpot/interface/BeamSpot.h"
88     #include "DataFormats/Common/interface/Handle.h"
89     #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
90     #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
91     #include "DataFormats/Candidate/interface/CompositePtrCandidate.h"
92     //trigger
93     #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuRegionalCand.h"
94     #include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTReadoutCollection.h"
95     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
96     #include "DataFormats/L1GlobalTrigger/interface/L1GtPsbWord.h"
97     #include "FWCore/Common/interface/TriggerNames.h"
98     #include "DataFormats/Common/interface/TriggerResults.h"
99     #include "DataFormats/HLTReco/interface/TriggerObject.h"
100     #include "DataFormats/HLTReco/interface/TriggerEvent.h"
101     #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
102     #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
103     //jpt
104     #include "DataFormats/JetReco/interface/JPTJetCollection.h"
105     #include "DataFormats/JetReco/interface/JPTJet.h"
106     //pf
107     #include "DataFormats/JetReco/interface/PFJetCollection.h"
108     #include "DataFormats/JetReco/interface/PFJet.h"
109     //vtx
110     // not sure if all that includes are needed .....
111     #include "DataFormats/TrackReco/interface/Track.h"
112     #include "DataFormats/TrackReco/interface/TrackFwd.h"
113     #include "DataFormats/VertexReco/interface/Vertex.h"
114     #include "DataFormats/VertexReco/interface/VertexFwd.h"
115     #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
116     #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
117     #include "DataFormats/Math/interface/deltaR.h"
118     #include "TTree.h"
119    
120 efe 1.2 //electron
121     //electron
122     #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
123     #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
124     #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
125     #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
126     #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
127     #include "DataFormats/EgammaReco/interface/BasicCluster.h"
128     #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
129     #include "DataFormats/PatCandidates/interface/Electron.h"
130    
131     //HFelectron
132     #include "DataFormats/EgammaReco/interface/HFEMClusterShape.h"
133     #include "DataFormats/EgammaReco/interface/HFEMClusterShapeAssociation.h"
134     #include "RecoEgamma/EgammaHFProducers/plugins/HFRecoEcalCandidateProducer.h"
135     #include "RecoEgamma/EgammaHFProducers/interface/HFRecoEcalCandidateAlgo.h"
136     #include "DataFormats/EgammaReco/interface/HFEMClusterShapeFwd.h"
137     #include "RecoEgamma/EgammaHFProducers/plugins/HFEMClusterProducer.h"
138     #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
139     #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
140     #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
141     #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
142     #include "DataFormats/EgammaReco/interface/SuperCluster.h"
143     #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
144    
145     #include "DataFormats/EgammaReco/interface/BasicCluster.h"
146     #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
147     #include <map>
148     #include "DataFormats/Common/interface/AssociationMap.h"
149    
150    
151 efe 1.1 class TH1F;
152     class TH2F;
153     class TStyle;
154     class TTree;
155    
156     //
157     // class decleration
158     //
159    
160     using namespace edm;
161     using namespace reco;
162     using namespace std;
163     using namespace ROOT::Math::VectorUtil;
164     using namespace HepMC;
165    
166    
167     class AFB : public edm::EDAnalyzer {
168     public:
169     explicit AFB(const edm::ParameterSet&);
170     ~AFB();
171    
172     const edm::ValueMap<double>& getValueMap(const edm::Event& iEvent,edm::InputTag& inputTag);
173    
174 efe 1.2 //void produce(const edm::Handle<SuperClusterCollection>& SuperClusters,
175     // const HFEMClusterShapeAssociationCollection& AssocShapes,
176     // RecoEcalCandidateCollection& RecoECand);
177    
178 efe 1.1
179     private:
180     // virtual void beginJob(const edm::EventSetup&) ;
181     virtual void beginJob();
182     virtual void analyze(const edm::Event&, const edm::EventSetup&);
183     virtual void endJob() ;
184    
185     bool IsMuMatchedToHLTMu ( const reco::Muon & , std::vector<reco::Particle> ,double ,double );
186    
187 efe 1.2
188    
189    
190    
191 efe 1.1 // ----------member data ---------------------------
192    
193     float DeltaPhi(float phi1, float phi2);
194     float DeltaR(float eta1, float eta2, float phi1, float phi2);
195     edm::InputTag muonTag_;
196     string CaloJetAlg;
197     edm::InputTag JPTAlg, JPTAlgL2L3;
198 efe 1.6 edm::InputTag trigTag_, triggerSummaryLabel_, hltTag_, hltTag2_, hltTag3_, hltTag4_, hltTag5_;
199 efe 1.1 string muonTrig_;
200 efe 1.6 string L3FilterName_, L3FilterName2_, L3FilterName3_, L3FilterName4_, L3FilterName5_;
201 efe 1.1 edm::Service<TFileService> fs;
202     TTree * myTree;
203     int event, run,lumi,bxnumber,realdata;
204     int hlt_trigger_fired;
205     int sort_index_for_mu_tree;
206     float RecMuonPt[50], RecMuonEta[50], RecMuonPhi[50],RecMuonPx[50], RecMuonPy[50], RecMuonPz[50], RecMuonE[50], RecMuonM[50], RecMuonGlobalType[50], RecMuonTrackerType[50], RecMuonStandAloneType[50], RecMuonIsoSumPt[50],RecMuonIsoRelative[50], RecMuonIsoCalComb[50], RecMuonglmuon_dxy[50], RecMuonglmuon_dz[50], RecMuonglmuon_normalizedChi2[50],RecMuonVx[50],RecMuonVy[50],RecMuonVz[50];
207 efe 1.6 int RecMuonglmuon_trackerHits[50],RecMuontkmuon_pixelhits[50],RecMuonglmuon_muonHits[50],RecNumberOfUsedStations[50],hltmatchedmuon[50],hltmatched_Dimuon[50];//,hltmatchedmuon2[50];
208 efe 1.2
209     float RecElec_Pt[50], RecElec_Px[50], RecElec_Py[50],RecElec_Pz[50],RecElec_eta[50],RecElec_phi[50],RecElec_GsfTrk_d0[50],RecElec_dr03TkSumPt[50],RecElec_dr03EcalRecHitSumEt[50],RecElec_dr03HcalTowerSumEt[50],RecElec_scSigmaIEtaIEta[50],RecElec_deltaPhiSuperClusterTrackAtVtx[50],RecElec_deltaEtaSuperClusterTrackAtVtx[50],RecElec_hadronicOverEm[50],RecElec_gsfTrack_numberOfLostHits[50];
210     int recelec_index,RecElec_Charge[50],RecElec_IsEB[50], RecElec_IsEE[50];
211    
212 efe 1.1 //particle information
213     int par_index, mom[50], daug[50];
214     float ParticlePt[50], ParticleEta[50], ParticlePhi[50], ParticlePx[50], ParticlePy[50], ParticlePz[50], ParticleE[50], ParticleM[50];
215     int ParticleId[50], ParticleStatus[50], ParticleMother[50][10], ParticleDaughter[50][10];
216     int id_tmp[20];
217     // int id_muon[20];
218     int RecMuonglmuon_charge[50];
219     //reco jets
220     int reco_jet;
221     //int jpt_c,cjpt_c
222     int pfNjets,pfNtracks;
223     float RecJetPt[50], RecJetEta[50], RecJetPhi[50], RecJetPx[50], RecJetPy[50], RecJetPz[50], RecJetE[50];
224     float RecCorrJetPt[50];
225     float JetEMF[50],JetN90[50],JetFHPD[50],JetFRBX[50];
226     float JPTPt[50], JPTEta[50], JPTPhi[50], JPTPx[50], JPTPy[50], JPTPz[50], JPTE[50];
227     float cJPTPt[50], cJPTEta[50], cJPTPhi[50], cJPTPx[50], cJPTPy[50], cJPTPz[50], cJPTE[50];
228     float PFjetEta[30], PFjetPhi[30],PFjetPt[30],PFCorrjetPt[30],PFjetCEMF[30],PFjetNEMF[30];
229     float PFjetTrkVZ[50][30],PFjetTrkPT[50][30];
230     float vtxZ[50],vtxZerr[50];
231     float vtxY[50],vtxYerr[50];
232     float vtxX[50],vtxXerr[50];
233     int vtxisValid[50],vtxisFake[50];
234     int nVertices,nGoodVertices;
235     int techTrigger[44];
236     //met
237     float caloMET, caloSET, pfMET, pfSET;
238     float caloMETX, pfMETX;
239     float caloMETY, pfMETY;
240     float muCorrMET, muCorrSET;
241    
242 efe 1.2
243    
244 efe 1.1 reco::helper::JetIDHelper *jetID;
245    
246     };
247    
248    
249    
250     AFB::AFB(const edm::ParameterSet& iConfig)
251     {
252     muonTag_ = iConfig.getParameter<edm::InputTag>("muonTag");
253     CaloJetAlg = iConfig.getParameter<string>("CaloJetAlg");
254     // JPTAlg = iConfig.getParameter<edm::InputTag>("JPTjets");
255     trigTag_ = iConfig.getParameter<edm::InputTag>("trigTag");
256     jetID = new reco::helper::JetIDHelper(iConfig.getParameter<ParameterSet>("JetIDParams"));
257     triggerSummaryLabel_ = iConfig.getParameter<edm::InputTag>("triggerSummaryLabel");
258     hltTag_ = iConfig.getParameter<edm::InputTag>("hltTag");
259 efe 1.4 hltTag2_ = iConfig.getParameter<edm::InputTag>("hltTag2");
260 efe 1.2 hltTag3_ = iConfig.getParameter<edm::InputTag>("hltTag3");
261 efe 1.6 hltTag4_ = iConfig.getParameter<edm::InputTag>("hltTag4");
262     hltTag5_ = iConfig.getParameter<edm::InputTag>("hltTag5");
263    
264 efe 1.1 L3FilterName_ = iConfig.getParameter<std::string>("L3FilterName");
265 efe 1.4 L3FilterName2_ = iConfig.getParameter<std::string>("L3FilterName2");
266 efe 1.2 L3FilterName3_ = iConfig.getParameter<std::string>("L3FilterName3");
267 efe 1.6 L3FilterName4_ = iConfig.getParameter<std::string>("L3FilterName4");
268     L3FilterName5_ = iConfig.getParameter<std::string>("L3FilterName5");
269 efe 1.2
270 efe 1.1 }
271    
272     AFB::~AFB()
273     {
274     }
275    
276    
277     // ------------ method called to for each event ------------
278     void
279     AFB::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
280     {
281    
282     bool singleTrigFlag1 = false;
283 efe 1.6 bool singleTrigFlag2 = false;
284     bool doubleTrigFlag1 = false;
285     bool doubleTrigFlag2 = false;
286 efe 1.1
287     Handle<GenParticleCollection> genParticles_h;
288     iEvent.getByLabel("genParticles", genParticles_h);
289     const GenParticleCollection* genParticles = genParticles_h.failedToGet () ? 0 : &*genParticles_h;
290    
291     Handle<MuonCollection> muons;
292     iEvent.getByLabel(muonTag_,muons);
293    
294     Handle<CaloJetCollection> caloJets;
295     iEvent.getByLabel(CaloJetAlg,caloJets);
296     CaloJetCollection::const_iterator c_jet;
297     // const JetCorrector* jCorrector = JetCorrector::getJetCorrector("L2L3JetCorrectorAK5Calo",iSetup);
298     const JetCorrector* jCorrector = JetCorrector::getJetCorrector("ak5CaloL2L3",iSetup);
299     const JetCorrector* PFjCorrector = JetCorrector::getJetCorrector("ak5PFL2L3",iSetup);
300    
301     /*
302     Handle<JPTJetCollection> jptJets;
303     iEvent.getByLabel(JPTAlg,jptJets);
304     */
305     /*
306     Handle<JPTJetCollection> corjptJets;
307     iEvent.getByLabel(JPTAlgL2L3,corjptJets);
308     */
309    
310     Handle< edm::View<reco::CaloMET> > caloMEThandle;
311     iEvent.getByLabel("met", caloMEThandle);
312    
313     Handle< edm::View<reco::PFMET> > pfMEThandle;
314     iEvent.getByLabel("pfMet", pfMEThandle);
315    
316     Handle< edm::View<reco::CaloMET> > muCorrMEThandle;
317     iEvent.getByLabel("corMetGlobalMuons", muCorrMEThandle);
318    
319     Handle<PFJetCollection> PFjets;
320     iEvent.getByLabel("ak5PFJets","",PFjets);
321     PFJetCollection::const_iterator i_jet;
322    
323 efe 1.2 //Electrons
324     edm::Handle<reco::GsfElectronCollection> elecs_h;
325     iEvent.getByLabel("gsfElectrons",elecs_h);
326     //iEvent.getByLabel("electronCollection_",elecs_h);
327     const GsfElectronCollection* elec = elecs_h.failedToGet() ? 0 : &*elecs_h;
328    
329     //HF electrons
330     edm::Handle<reco::RecoEcalCandidateCollection> hfelecs_h;
331     iEvent.getByLabel("hfRecoEcalCandidate",hfelecs_h);
332     const RecoEcalCandidateCollection* hfelec = hfelecs_h.failedToGet() ? 0 : &*hfelecs_h;
333    
334     Handle<SuperClusterCollection> SCH_h;
335     iEvent.getByLabel("hfEMClusters",SCH_h);
336     const SuperClusterCollection* SCH = SCH_h.failedToGet() ? 0 : &*SCH_h; //???
337     Handle<HFEMClusterShapeAssociationCollection> ASH_h;
338     iEvent.getByLabel("hfEMClusters",ASH_h);
339     const HFEMClusterShapeAssociationCollection* ASH = ASH_h.failedToGet() ? 0 : &*ASH_h;//??
340    
341    
342 efe 1.1
343     Handle<BeamSpot> beamSpotHandle;
344     if (!iEvent.getByLabel(InputTag("offlineBeamSpot"), beamSpotHandle)) {
345     cout<< ">>> No beam spot found !!!" <<endl;
346     return;
347     }
348    
349     Handle<VertexCollection> pvHandle;
350     iEvent.getByLabel("offlinePrimaryVertices", pvHandle);
351    
352     nVertices = pvHandle->size();
353    
354     const VertexCollection vertexColl = *(pvHandle.product());
355    
356     nGoodVertices = 0;
357     for (VertexCollection::const_iterator vtx = vertexColl.begin(); vtx!= vertexColl.end(); ++vtx){
358     if (vtx->isValid() && !vtx->isFake()) ++nGoodVertices;
359     }
360    
361     // if (nVertices != nGoodVertices) cout<<"@@@@@@@@@@@@@ ------> "<<nVertices<<" "<<nGoodVertices<<endl;
362    
363     reco::VertexCollection vtxs = *pvHandle;
364     for (int iv = 0;iv<nVertices;iv++){
365     vtxX[iv] = vtxs[iv].x();
366     vtxY[iv] = vtxs[iv].y();
367     vtxZ[iv] = vtxs[iv].z();
368     vtxXerr[iv] = vtxs[iv].xError();
369     vtxYerr[iv] = vtxs[iv].yError();
370     vtxZerr[iv] = vtxs[iv].zError();
371     vtxisValid[iv] = vtxs[iv].isValid();
372     vtxisFake[iv] = vtxs[iv].isFake();
373     }
374    
375     Handle< L1GlobalTriggerReadoutRecord > gtRecord;
376     iEvent.getByLabel( "gtDigis", gtRecord);
377     const TechnicalTriggerWord tWord = gtRecord->technicalTriggerWord();
378    
379     Handle<TriggerResults> triggerResults;
380     if (!iEvent.getByLabel(trigTag_, triggerResults)){
381     cout<<" Hoooop"<<endl;
382     LogError("") << ">>> TRIGGER collection does not exist !!!";
383     };
384     // const edm::TriggerNames & triggerNames = iEvent.triggerNames(*triggerResults);
385     // int hlt_trigger_fired = 0;
386     // int itrig1 = triggerNames.triggerIndex(muonTrig_);
387     // if (triggerResults->accept(itrig1)) hlt_trigger_fired = 1;
388    
389 efe 1.3
390     for (int ee = 0;ee<44;ee++){
391     techTrigger[ee] = tWord.at(ee);
392     }
393    
394     event = iEvent.id().event();
395     run = iEvent.id().run();
396     lumi = iEvent.luminosityBlock();
397     bxnumber = iEvent.bunchCrossing();
398     realdata = iEvent.isRealData();
399    
400 efe 1.1 //---------------------For trigger matching-----------------------------------
401    
402     Handle<trigger::TriggerEvent> triggerObj;
403 efe 1.3 iEvent.getByLabel(triggerSummaryLabel_,triggerObj);
404 efe 1.1 const trigger::TriggerObjectCollection & toc(triggerObj->getObjects());
405     size_t nMuHLT =0;
406     std::vector<reco::Particle> HLTMuMatched;
407 efe 1.6 //for dimuons
408     size_t nDiMuHLT =0;
409     std::vector<reco::Particle> HLTDiMuMatched;
410     //--end of for dimuons
411 efe 1.1 for ( size_t ia = 0; ia < triggerObj->sizeFilters(); ++ ia) {
412     std::string fullname = triggerObj->filterTag(ia).encode();
413 efe 1.3 // cout<<"full name: "<<fullname<<endl;
414 efe 1.1 std::string name;
415     size_t p = fullname.find_first_of(':');
416     if ( p != std::string::npos) {
417     name = fullname.substr(0, p);
418     }
419     else {
420     name = fullname;
421     }
422     if ( &toc !=0 ) {
423     const trigger::Keys & k = triggerObj->filterKeys(ia);
424     for (trigger::Keys::const_iterator ki = k.begin(); ki !=k.end(); ++ki ) {
425 efe 1.2 // if (name == L3FilterName_ ) {
426 efe 1.4 if ( (run <= 146111 && name == L3FilterName_) || (run >=147196 && run <148108 && name == L3FilterName2_) || ( run >=148108 && name == L3FilterName3_) ) {
427     cout<<name << "=?" << L3FilterName_<<" "<<L3FilterName2_<<" "<<L3FilterName3_<<endl;
428 efe 1.1 HLTMuMatched.push_back(toc[*ki].particle());
429     nMuHLT++;
430     }
431 efe 1.6 if ( (run < 147196 && name == L3FilterName4_) || (run >= 147196 && name == L3FilterName5_)){
432     cout<<name << "=?" << L3FilterName4_<<" "<<L3FilterName5_<<endl;
433     HLTDiMuMatched.push_back(toc[*ki].particle());
434     nDiMuHLT++;
435     }
436 efe 1.1 }
437     }
438     }
439    
440     //----------------------------------------------------------------------------
441    
442 efe 1.3
443 efe 1.1 //--------------------met-------------------------------------------------
444     caloMET = (caloMEThandle->front()).et();
445     caloSET = (caloMEThandle->front()).sumEt();
446     pfMET = (pfMEThandle->front()).et();
447     pfSET = (pfMEThandle->front()).sumEt();
448    
449     caloMETX = (caloMEThandle->front()).px();
450     caloMETY = (caloMEThandle->front()).py();
451    
452     pfMETX = (pfMEThandle->front()).px();
453     pfMETY = (pfMEThandle->front()).py();
454    
455     muCorrMET = (muCorrMEThandle->front()).et();
456     muCorrSET = (muCorrMEThandle->front()).sumEt();
457     //------------------------------------------------------------------------
458    
459     par_index = 0;
460     for (int i=0;i<50;i++){
461     ParticlePt[i] = -999;
462     ParticleEta[i] = -999;
463     ParticlePhi[i] = -999;
464     ParticlePx[i] = -999;
465     ParticlePy[i] = -999;
466     ParticlePz[i] = -999;
467     ParticleE[i] = -999;
468     ParticleM[i] = -999;
469     ParticleId[i] = -999;
470     ParticleStatus[i] = -999;
471     for (int j=0;j<10;j++){
472     ParticleMother[i][j] = -999;
473     ParticleDaughter[i][j] = -999;
474     }
475     }
476    
477     //--------------------particles-------------------------------------------
478    
479     if (!realdata && genParticles){
480     par_index = 0;
481     for (size_t i=0; i<genParticles->size(); ++i){
482     //const GenParticle & p = (*genParticles)[i];
483     const Candidate & p = (*genParticles)[i];
484     int id = p.pdgId();
485     int st = p.status();
486     if (st!=3) continue;
487     ParticlePt[par_index] = p.pt();
488     ParticleEta[par_index] = p.eta();
489     ParticlePhi[par_index] = p.phi();
490     ParticlePx[par_index] = p.px();
491     ParticlePy[par_index] = p.py();
492     ParticlePz[par_index] = p.pz();
493     ParticleE[par_index] =p.energy();
494     ParticleM[par_index] =p.mass();
495     ParticleId[par_index] = id;
496     ParticleStatus[par_index] = st;
497     //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
498    
499     //cout<<i+1<<" "<<st<<" "<<id<<" "<<endl;
500     mom[par_index] = p.numberOfMothers();
501     for (uint d=0; d<p.numberOfMothers(); d++) {
502     // cout << " mother id: " << p.mother(d)->pdgId() << " ";
503     // mom = d;
504     ParticleMother[par_index][d] = p.mother(d)->pdgId();
505     }
506     daug[par_index] = p.numberOfDaughters();
507     for (uint d=0; d<p.numberOfDaughters(); d++) {
508     if (p.daughter(d)->status() == 3){
509     // cout << " daughter id: " << p.daughter(d)->pdgId() << " ";
510     // daug = d;
511     ParticleDaughter[par_index][d] = p.daughter(d)->pdgId();
512     }
513     }
514     // cout<<" "<<endl;
515     //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
516     ++par_index;
517     // cout<< p.pdgId() << " "<< p.mass() << ", "<< p.status() << endl;
518     // cout<<i<<" "<<st<<" "<<id<<" "<<p.eta()<<" "<<p.mass()<<endl;
519     }
520     }
521    
522     //------------------------------------------------------------------------
523    
524     //--------------------------MUONS---------------------------------------------------------------------------------------------------------------------
525    
526     int reco_muon = 0;
527     float rec_eta_muon[50] = {};
528     float rec_phi_muon[50] = {};
529     float rec_pt_muon[50] = {};
530     float rec_px_muon[50] = {};
531     float rec_py_muon[50] = {};
532     float rec_pz_muon[50] = {};
533     float rec_vx_muon[50] = {};
534     float rec_vy_muon[50] = {};
535     float rec_vz_muon[50] = {};
536     float rec_e_muon[50] = {};
537     float rec_m_muon[50] = {};
538     int rec_n_used_sta[50] = {};
539     float rec_typeGlobal_muon[50] = {};
540     float rec_typeTracker_muon[50] = {};
541     float rec_typeStandAlone_muon[50] = {};
542     float rec_iso_sumpt[50] = {};
543     float rec_iso_relative[50] = {};
544     float rec_iso_CalComb[50] = {};
545     float glmuon_dxy[50] = {};
546     float glmuon_dz[50] = {};
547     float glmuon_normalizedChi2[50] = {};
548     int glmuon_trackerHits[50] = {};
549     int tkmuon_pixelhits[50] = {};
550     int glmuon_muonHits[50] = {};
551     int glmuon_charge[50] = {};
552     int hlt_tmp[50] = {};
553 efe 1.6 // int hlt_tmp2[50] = {};
554     int hlt_tmp_dimuon[50] = {};
555     // int hlt_tmp_dimuon2[50] = {};
556 efe 1.5 for (int gg = 0;gg<50;gg++){
557     rec_eta_muon[gg] = 0;
558     rec_phi_muon[gg] = 0;
559     rec_pt_muon[gg] = 0;
560     rec_px_muon[gg] = 0;
561     rec_py_muon[gg] = 0;
562     rec_pz_muon[gg] = 0;
563     rec_vx_muon[gg] = 0;
564     rec_vy_muon[gg] = 0;
565     rec_vz_muon[gg] = 0;
566     rec_e_muon[gg] = 0;
567     rec_m_muon[gg] = 0;
568     rec_n_used_sta[gg] = 0;
569     rec_typeGlobal_muon[gg] = 0;
570     rec_typeTracker_muon[gg] = 0;
571     rec_typeStandAlone_muon[gg] = 0;
572     rec_iso_sumpt[gg] = 0;
573     rec_iso_relative[gg] = 0;
574     rec_iso_CalComb[gg] = 0;
575     glmuon_dxy[gg] = 0;
576     glmuon_dz[gg] = 0;
577     glmuon_normalizedChi2[gg] = 0;
578     glmuon_trackerHits[gg] = 0;
579     tkmuon_pixelhits[gg] = 0;
580     glmuon_muonHits[gg] = 0;
581     glmuon_charge[gg] = 0;
582 efe 1.6 hlt_tmp[gg] = 0;
583     // hlt_tmp2[gg] = 0;
584 efe 1.5 }
585 efe 1.1 MuonCollection::const_iterator muon;
586     std::vector<reco::Muon> highPtGlbMuons;
587     for (unsigned int i=0; i<muons->size(); i++ ){
588     const reco::Muon & mu = muons->at(i);
589     // if (mu.isGlobalMuon()) highPtGlbMuons.push_back(mu);
590     if (mu.isGlobalMuon() && mu.isTrackerMuon()){
591     highPtGlbMuons.push_back(mu);
592     rec_typeGlobal_muon[reco_muon] = mu.isGlobalMuon();
593     rec_typeTracker_muon[reco_muon] = mu.isTrackerMuon();
594     if (mu.isStandAloneMuon()) rec_typeStandAlone_muon[reco_muon] = mu.isStandAloneMuon();
595     reco::TrackRef glmuon = mu.globalTrack();//OK
596     reco::TrackRef tkmuon = mu.innerTrack();//OK
597     glmuon_dxy[reco_muon] = glmuon->dxy(beamSpotHandle->position());//OK
598     glmuon_dz[reco_muon] = glmuon->dz(beamSpotHandle->position());//probably OK but not used anyway.
599     glmuon_normalizedChi2[reco_muon] = glmuon->normalizedChi2();//OK
600     glmuon_trackerHits[reco_muon] = tkmuon->hitPattern().numberOfValidTrackerHits();//OK
601     tkmuon_pixelhits[reco_muon] = tkmuon->hitPattern().numberOfValidPixelHits();//OK
602     glmuon_muonHits[reco_muon] = glmuon->hitPattern().numberOfValidMuonHits();//OK
603     glmuon_charge[reco_muon] = glmuon->charge();
604     rec_eta_muon[reco_muon] = mu.eta();
605     rec_phi_muon[reco_muon] = mu.phi();
606     rec_pt_muon[reco_muon] = mu.pt();
607     rec_px_muon[reco_muon] = mu.px();
608     rec_py_muon[reco_muon] = mu.py();
609     rec_pz_muon[reco_muon] = mu.pz();
610     rec_vx_muon[reco_muon] = mu.vx();
611     rec_vy_muon[reco_muon] = mu.vy();
612     rec_vz_muon[reco_muon] = mu.vz();
613     rec_e_muon[reco_muon] = mu.energy();
614     rec_m_muon[reco_muon] = mu.mass();
615     rec_n_used_sta[reco_muon] = mu.numberOfMatches();
616     // id_tmp[reco_muon] = muon::isGoodMuon(*muon,muon::GlobalMuonPromptTight);//id
617     rec_iso_sumpt[reco_muon] = mu.isolationR03().sumPt;
618     rec_iso_relative[reco_muon] = rec_iso_sumpt[reco_muon]/mu.pt();//OK
619     rec_iso_CalComb[reco_muon] = mu.isolationR03().emEt + mu.isolationR03().hadEt;
620     /*
621     reco::Muon muon1 = highPtGlbMuons[i];//PROBLEM IS HERE!!!!!!!!!!!!!!!!!!!!!!!
622     math::XYZTLorentzVector mu1(muon1.p4());
623     singleTrigFlag1 = IsMuMatchedToHLTMu ( muon1, HLTMuMatched , 0.2, 1.0 );
624     hlt_tmp[reco_muon] = int(singleTrigFlag1);
625     */
626     ++reco_muon;
627     }
628     }
629    
630     //==============================================================
631     unsigned int nHighPtGlbMu = highPtGlbMuons.size();
632     if (nHighPtGlbMu > 0){
633 efe 1.6 /*
634 efe 1.1 for(unsigned int i =0 ; i < nHighPtGlbMu ; i++) {
635     reco::Muon muon1 = highPtGlbMuons[i];
636     singleTrigFlag1 = IsMuMatchedToHLTMu ( muon1, HLTMuMatched , 0.2, 1.0 );
637     hlt_tmp[i] = int(singleTrigFlag1);
638     }
639 efe 1.6 */
640     if (nHighPtGlbMu>1 ){
641     for(unsigned int i =0 ; i < nHighPtGlbMu ; i++) {
642     reco::Muon muon1 = highPtGlbMuons[i];
643     for (unsigned int j =i+1; j <nHighPtGlbMu ; ++j ){
644     reco::Muon muon2 = highPtGlbMuons[j];
645     if (muon1.charge() == muon2.charge()) continue;
646     if (nMuHLT){
647     singleTrigFlag1 = IsMuMatchedToHLTMu ( muon1, HLTMuMatched , 0.2, 1.0 );
648     singleTrigFlag2 = IsMuMatchedToHLTMu ( muon2, HLTMuMatched , 0.2, 1.0 );
649     hlt_tmp[i] = int(singleTrigFlag1);
650     hlt_tmp[j] = int(singleTrigFlag2);
651     // cout<<"SINGLE i,j hlt tmp1,2 : "<< i<<" "<<j<<" "<<hlt_tmp[i] <<" "<< hlt_tmp[j]<<endl;
652     }
653     if (nDiMuHLT){
654     doubleTrigFlag1 = IsMuMatchedToHLTMu ( muon1, HLTDiMuMatched , 0.2, 1.0 );
655     doubleTrigFlag2 = IsMuMatchedToHLTMu ( muon2, HLTDiMuMatched , 0.2, 1.0 );
656     hlt_tmp_dimuon[i] = int(doubleTrigFlag1);
657     hlt_tmp_dimuon[j] = int(doubleTrigFlag2);
658     // cout<<"DOUBLE i,j hlt tmp dimuon 1,2 : "<< i<<" "<<j<<" "<<hlt_tmp_dimuon1[i] <<" "<< hlt_tmp_dimuon2[j]<<endl;
659     }
660     }
661     }
662     }
663     for(unsigned int i =0 ; i < nHighPtGlbMu ; i++){
664     if (nMuHLT) cout<<"SINGLE hlt tmp: "<< i<<" "<<hlt_tmp[i] <<endl;
665     }
666     for(unsigned int i =0 ; i < nHighPtGlbMu ; i++){
667     if (nDiMuHLT) cout<<"DOUBLE hlt tmp: "<< i<<" "<<hlt_tmp_dimuon[i] <<endl;
668     }
669 efe 1.1 }
670     //==============================================================
671    
672     // cout<<"Number of Muons: "<<muons->size()<<endl;
673     /*
674     //----old one---------
675     for(muon = muons->begin(); muon != muons->end(); ++muon){
676     if (muon->isGlobalMuon()) rec_typeGlobal_muon[reco_muon] = muon->isGlobalMuon();
677     if (muon->isTrackerMuon()) rec_typeTracker_muon[reco_muon] = muon->isTrackerMuon();
678     if (muon->isStandAloneMuon()) rec_typeStandAlone_muon[reco_muon] = muon->isStandAloneMuon();
679     if (muon->isGlobalMuon() && muon->isTrackerMuon()){
680     reco::TrackRef glmuon = muon->globalTrack();//OK
681     reco::TrackRef tkmuon = muon->innerTrack();//OK
682     glmuon_dxy[reco_muon] = glmuon->dxy(beamSpotHandle->position());//OK
683     glmuon_dz[reco_muon] = glmuon->dz(beamSpotHandle->position());//probably OK but not used anyway.
684     glmuon_normalizedChi2[reco_muon] = glmuon->normalizedChi2();//OK
685     glmuon_trackerHits[reco_muon] = tkmuon->hitPattern().numberOfValidTrackerHits();//OK
686     tkmuon_pixelhits[reco_muon] = tkmuon->hitPattern().numberOfValidPixelHits();//OK
687     glmuon_muonHits[reco_muon] = glmuon->hitPattern().numberOfValidMuonHits();//OK
688     glmuon_charge[reco_muon] = glmuon->charge();
689     rec_eta_muon[reco_muon] = muon->eta();
690     rec_phi_muon[reco_muon] = muon->phi();
691     rec_pt_muon[reco_muon] = muon->pt();
692     rec_px_muon[reco_muon] = muon->px();
693     rec_py_muon[reco_muon] = muon->py();
694     rec_pz_muon[reco_muon] = muon->pz();
695     rec_e_muon[reco_muon] = muon->energy();
696     rec_m_muon[reco_muon] = muon->mass();
697     rec_n_used_sta[reco_muon] = muon->numberOfMatches();
698     id_tmp[reco_muon] = muon::isGoodMuon(*muon,muon::GlobalMuonPromptTight);//id
699     rec_iso_sumpt[reco_muon] = muon->isolationR03().sumPt;
700     rec_iso_relative[reco_muon] = rec_iso_sumpt[reco_muon]/muon->pt();//OK
701     rec_iso_CalComb[reco_muon] = muon->isolationR03().emEt + muon->isolationR03().hadEt;
702     ++reco_muon;
703     }
704     }
705     */
706     //-------------
707    
708     for (int mm=0;mm<reco_muon;++mm){
709     RecMuonPt[mm] = 0;
710     RecMuonEta[mm] = 0;
711     RecMuonPhi[mm] = 0;
712     RecMuonPx[mm] = 0;
713     RecMuonPy[mm] = 0;
714     RecMuonPz[mm] = 0;
715     RecMuonVx[mm] = 0;
716     RecMuonVy[mm] = 0;
717     RecMuonVz[mm] = 0;
718     RecMuonE[mm] = 0;
719     RecMuonM[mm] = 0;
720     RecMuonGlobalType[mm] = 0;
721     RecMuonStandAloneType[mm] = 0;
722     RecMuonTrackerType[mm] = 0;
723     RecMuonIsoSumPt[mm] = 0;
724     RecNumberOfUsedStations[mm] = 0;
725     RecMuonIsoRelative[mm] = 0;
726     RecMuonIsoCalComb[mm] = 0;
727     RecMuonglmuon_dxy[mm] = 0;
728     RecMuonglmuon_dz[mm] = 0;
729     RecMuonglmuon_normalizedChi2[mm] = 0;
730     RecMuonglmuon_trackerHits[mm] = 0;
731     RecMuontkmuon_pixelhits[mm] = 0;
732     RecMuonglmuon_muonHits[mm] = 0;
733     RecMuonglmuon_charge[mm] = 0;
734     // id_muon[mm] = 0;
735     }
736    
737     int sorted_muon_index[10]={};
738     int sort_reco_muon = 0;
739     sort_index_for_mu_tree = 0;
740     TMath::Sort(reco_muon,rec_pt_muon,sorted_muon_index);
741     for (int i=0;i<reco_muon;++i){
742     sort_reco_muon = sorted_muon_index[sort_index_for_mu_tree];
743     RecMuonPt[sort_index_for_mu_tree] = rec_pt_muon[sort_reco_muon];
744     RecMuonEta[sort_index_for_mu_tree] = rec_eta_muon[sort_reco_muon];
745     RecMuonPhi[sort_index_for_mu_tree] = rec_phi_muon[sort_reco_muon];
746     RecMuonPx[sort_index_for_mu_tree] = rec_px_muon[sort_reco_muon];
747     RecMuonPy[sort_index_for_mu_tree] = rec_py_muon[sort_reco_muon];
748     RecMuonPz[sort_index_for_mu_tree] = rec_pz_muon[sort_reco_muon];
749     RecMuonVx[sort_index_for_mu_tree] = rec_vx_muon[sort_reco_muon];
750     RecMuonVy[sort_index_for_mu_tree] = rec_vy_muon[sort_reco_muon];
751     RecMuonVz[sort_index_for_mu_tree] = rec_vz_muon[sort_reco_muon];
752     RecMuonE[sort_index_for_mu_tree] = rec_e_muon[sort_reco_muon];
753     RecMuonM[sort_index_for_mu_tree] = rec_m_muon[sort_reco_muon];
754     RecMuonGlobalType[sort_index_for_mu_tree] = rec_typeGlobal_muon[sort_reco_muon];
755     RecMuonTrackerType[sort_index_for_mu_tree] = rec_typeTracker_muon[sort_reco_muon];
756     RecMuonStandAloneType[sort_index_for_mu_tree] = rec_typeStandAlone_muon[sort_reco_muon];
757     RecMuonIsoSumPt[sort_index_for_mu_tree] = rec_iso_sumpt[sort_reco_muon];
758     RecNumberOfUsedStations[sort_index_for_mu_tree] = rec_n_used_sta[sort_reco_muon];
759     RecMuonIsoRelative[sort_index_for_mu_tree] = rec_iso_relative[reco_muon];
760     RecMuonIsoCalComb[sort_index_for_mu_tree] = rec_iso_CalComb[sort_reco_muon];
761     RecMuonglmuon_dxy[sort_index_for_mu_tree] = glmuon_dxy[sort_reco_muon];
762     RecMuonglmuon_dz[sort_index_for_mu_tree] = glmuon_dz[sort_reco_muon];
763     RecMuonglmuon_normalizedChi2[sort_index_for_mu_tree] = glmuon_normalizedChi2[sort_reco_muon];
764     RecMuonglmuon_trackerHits[sort_index_for_mu_tree] = glmuon_trackerHits[sort_reco_muon];
765     RecMuontkmuon_pixelhits[sort_index_for_mu_tree] = tkmuon_pixelhits[sort_reco_muon];
766     RecMuonglmuon_muonHits[sort_index_for_mu_tree] = glmuon_muonHits[sort_reco_muon];
767     RecMuonglmuon_charge[sort_index_for_mu_tree] = glmuon_charge[sort_reco_muon];
768     hltmatchedmuon[sort_index_for_mu_tree] = hlt_tmp[sort_reco_muon];
769 efe 1.6 hltmatched_Dimuon[sort_index_for_mu_tree] = hlt_tmp_dimuon[sort_reco_muon];
770     // hltmatchedmuon2[sort_index_for_mu_tree] = hlt_tmp2[sort_reco_muon];
771     cout<<"INDEX AND MATCHED MUONS IN THE TREE and pt(SINGLE): "<<sort_index_for_mu_tree<<" "<<hltmatchedmuon[sort_index_for_mu_tree]<<" "<<RecMuonPt[sort_index_for_mu_tree]<<endl;
772     cout<<"INDEX AND MATCHED MUONS IN THE TREE and pt (DOUBLE): "<<sort_index_for_mu_tree<<" "<<hltmatched_Dimuon[sort_index_for_mu_tree]<<" "<<RecMuonPt[sort_index_for_mu_tree]<<endl;
773    
774 efe 1.1 /*
775     cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
776     cout<<"Written pt: "<<RecMuonPt[sort_index_for_mu_tree]<<endl;
777     cout<<"matched ?: "<<hltmatchedmuon[sort_index_for_mu_tree]<<endl;
778     cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
779     */
780     // id_muon[sort_index_for_mu_tree] = id_tmp[sort_reco_muon];
781     ++sort_index_for_mu_tree;
782     }
783    
784     //-----------------------END MUONS--------------------------------------------------------------------------------------------------------------------
785    
786 efe 1.2 //--------------------------ELECTRONS-----------------------------------------------------------------------------------------------------------------
787     for (int jj = 0; jj<50; ++jj){
788     RecElec_Pt[jj] = -99.;
789     RecElec_Px[jj] = -99.;
790     RecElec_Py[jj] = -99.;
791     RecElec_Pz[jj] = -99.;
792     RecElec_eta[jj] = -99.;
793     RecElec_phi[jj] = -99.;
794     RecElec_Charge[jj] = -99.;
795     RecElec_GsfTrk_d0[jj] = -99.;
796     RecElec_IsEB[jj] = -99.;
797     RecElec_IsEE[jj] = -99.;
798     RecElec_dr03TkSumPt[jj] = -99.;
799     RecElec_dr03EcalRecHitSumEt[jj] = -99.;
800     RecElec_dr03HcalTowerSumEt[jj] = -99.;
801     RecElec_scSigmaIEtaIEta[jj] = -99.;
802     RecElec_deltaPhiSuperClusterTrackAtVtx[jj] = -99.;
803     RecElec_deltaEtaSuperClusterTrackAtVtx[jj] = -99.;
804     RecElec_hadronicOverEm[jj] = -99.;
805     RecElec_gsfTrack_numberOfLostHits[jj] = -99.;
806    
807     }
808     recelec_index = 0;
809     for(GsfElectronCollection::const_iterator el = elec->begin(); el != elec->end(); ++el) {
810     RecElec_Pt[recelec_index] = el->pt();
811     RecElec_Px[recelec_index] = el->px();
812     RecElec_Py[recelec_index] = el->py();
813     RecElec_Pz[recelec_index] = el->pz();
814     RecElec_eta[recelec_index] = el->eta();
815     RecElec_phi[recelec_index] = el->phi();
816     RecElec_Charge[recelec_index] = el->charge();
817     RecElec_GsfTrk_d0[recelec_index] = el->gsfTrack()->d0();
818     RecElec_IsEB[recelec_index] = el->isEB();
819     RecElec_IsEE[recelec_index] = el->isEE();
820     RecElec_dr03TkSumPt[recelec_index] = el->dr03TkSumPt();
821     RecElec_dr03EcalRecHitSumEt[recelec_index] = el->dr03EcalRecHitSumEt();
822     RecElec_dr03HcalTowerSumEt[recelec_index] = el->dr03HcalTowerSumEt();
823     RecElec_scSigmaIEtaIEta[recelec_index] = el->scSigmaIEtaIEta();
824     RecElec_deltaPhiSuperClusterTrackAtVtx[recelec_index] = el->deltaPhiSuperClusterTrackAtVtx();
825     RecElec_deltaEtaSuperClusterTrackAtVtx[recelec_index] = el->deltaEtaSuperClusterTrackAtVtx();
826     RecElec_hadronicOverEm[recelec_index] = el->hadronicOverEm();
827     RecElec_gsfTrack_numberOfLostHits[recelec_index] = el->gsfTrack()->numberOfLostHits();
828     ++recelec_index;
829     }
830    
831     for(reco::RecoEcalCandidateCollection::const_iterator hfe = hfelec->begin(); hfe != hfelec->end(); ++hfe){
832     // cout<<hfe->et()<<" "<<hfe->px()<<" "<<hfe->py()<<" "<<hfe->pz()<<" "<<hfe->energy()<<endl;
833     }
834    
835    
836     //--------------------------END ELECTRONS-----------------------------------------------------------------------------------------------------------------
837    
838    
839 efe 1.1
840     //-----------------------RECO JETS--------------------------------------------------------------------------------------------------------------------
841     reco_jet = 0;
842     for( CaloJetCollection::const_iterator jet = caloJets->begin(); jet != caloJets->end(); ++ jet ) {
843     const math::XYZTLorentzVector theJet = jet->p4();
844     RecJetPt[reco_jet] = jet->pt();
845     RecCorrJetPt[reco_jet] = jCorrector->correction(theJet)*jet->pt();
846     RecJetEta[reco_jet] = jet->eta();
847     RecJetPhi[reco_jet] = jet->phi();
848     RecJetPx[reco_jet] = jet->px();
849     RecJetPy[reco_jet] = jet->py();
850     RecJetPz[reco_jet] = jet->pz();
851     RecJetE[reco_jet] = jet->energy();
852     JetEMF[reco_jet] = jet->emEnergyFraction();
853    
854     jetID->calculate(iEvent,*jet);
855     JetN90[reco_jet] = jetID->n90Hits();
856     JetFHPD[reco_jet] = jetID->fHPD();
857     JetFRBX[reco_jet] = jetID->fRBX();
858    
859     ++reco_jet;
860     }
861    
862     //-----------------------END RECOJETS-----------------------------------------------------------------------------------------------------------------
863    
864     //----------------------JPT JETS-----------------------------------------------------------------------------------------------------------------
865     /*
866     jpt_c = 0;
867     for (int i=0;i<50;i++){
868     JPTPt[i] = -9999;
869     JPTEta[i] = -9999;
870     JPTPhi[i] = -9999;
871     JPTPx[i] = -9999;
872     JPTPy[i] = -9999;
873     JPTPz[i] = -9999;
874     JPTE[i] = -9999;
875     cJPTPt[i] = -9999;
876     cJPTEta[i] = -9999;
877     cJPTPhi[i] = -9999;
878     cJPTPx[i] = -9999;
879     cJPTPy[i] = -9999;
880     cJPTPz[i] = -9999;
881     cJPTE[i] = -9999;
882     }
883     for( JPTJetCollection::const_iterator jjet = jptJets->begin(); jjet != jptJets->end(); ++ jjet ) {
884     JPTPt[jpt_c] = jjet->pt();
885     JPTEta[jpt_c] = jjet->eta();
886     JPTPhi[jpt_c] = jjet->phi();
887     JPTPx[jpt_c] = jjet->px();
888     JPTPy[jpt_c] = jjet->py();
889     JPTPz[jpt_c] = jjet->pz();
890     JPTE[jpt_c] = jjet->energy();
891     ++jpt_c;
892     }
893     */
894     /*
895     cjpt_c = 0;
896     for( CaloJetCollection::const_iterator cjjet = corjptJets->begin(); cjjet != corjptJets->end(); ++ cjjet ) {//L2L3 corr JPT
897     cJPTPt[cjpt_c] = cjjet->pt();
898     cJPTEta[cjpt_c] = cjjet->eta();
899     cJPTPhi[cjpt_c] = cjjet->phi();
900     cJPTPx[cjpt_c] = cjjet->px();
901     cJPTPy[cjpt_c] = cjjet->py();
902     cJPTPz[cjpt_c] = cjjet->pz();
903     cJPTE[cjpt_c] = cjjet->energy();
904     ++cjpt_c;
905     }
906     */
907    
908    
909     //----------------------END JPT JETS-----------------------------------------------------------------------------------------------------------------
910    
911    
912    
913     int NJets = 20;
914     //----------------------------PF jets-------------------------------------------------------------------------------------
915     for (int kkk = 0;kkk<50;kkk++){
916     for (int jjj = 0;jjj<30;jjj++){
917     PFjetTrkVZ[kkk][jjj] = -99;
918     PFjetTrkPT[kkk][jjj] = -99;
919     }
920     }
921     pfNjets=0;
922     for(i_jet = PFjets->begin(); i_jet != PFjets->end() && pfNjets < NJets; i_jet++){
923     const math::XYZTLorentzVector theJet = i_jet->p4();
924     PFjetEta[pfNjets]=i_jet->eta();
925     PFjetPhi[pfNjets]=i_jet->phi();
926     PFjetPt[pfNjets]=i_jet->pt();
927     PFCorrjetPt[pfNjets] = PFjCorrector->correction(theJet)*i_jet->pt();
928     PFjetCEMF[pfNjets]=i_jet->chargedEmEnergyFraction();
929     PFjetNEMF[pfNjets]=i_jet->neutralEmEnergyFraction();
930     // cout<<PFjetEta[pfNjets]<<" "<<PFjetPt[pfNjets]<<" "<<PFCorrjetPt[pfNjets]<<endl;
931     // ---- This accesses to the vertex Z of the track-base constituents of pfjets
932     pfNtracks=0;
933     const reco::TrackRefVector &tracks = i_jet->getTrackRefs();
934     for (reco::TrackRefVector::const_iterator iTrack = tracks.begin();iTrack != tracks.end(); ++iTrack) {
935     PFjetTrkVZ[pfNjets][pfNtracks] = (**iTrack).vz();
936     PFjetTrkPT[pfNjets][pfNtracks] = (**iTrack).pt();
937     pfNtracks++;
938     }
939     pfNjets++;
940     }
941     //-----------------------------end of pf jets--------------------------------------------------------------------------
942    
943 efe 1.2 // cout<<"clusShape.e9e25() "<<clusShape.e9e25()<<" "<<"clusShape.eCOREe9() "<<clusShape.eCOREe9()<<" "<<"(clusShape.eSeL() "<<(clusShape.eSeL()<<endl;
944 efe 1.1
945     myTree->Fill();//!!!!!!
946    
947     }
948    
949    
950     // ------------ method called once each job just before starting event loop ------------
951     //void AFB::beginJob(const edm::EventSetup&)
952     void AFB::beginJob()
953     {
954     TFileDirectory TestDir = fs->mkdir("test");
955     myTree = new TTree("MuonTree","MuonTree");
956     myTree->Branch("event", &event, "event/I");
957     myTree->Branch("run", &run, "run/I");
958     myTree->Branch("lumi", &lumi, "lumi/I");
959     myTree->Branch("bxnumber", &bxnumber, "bxnumber/I");
960     myTree->Branch("realdata", &realdata, "realdata/I");
961    
962     //muon
963     myTree->Branch("hlt_trigger_fired",&hlt_trigger_fired,"hlt_trigger_fired/I");
964     myTree->Branch("sort_index_for_mu_tree",&sort_index_for_mu_tree,"sort_index_for_mu_tree/I");
965     myTree->Branch("RecMuonPt",RecMuonPt,"RecMuonPt[sort_index_for_mu_tree]/F");
966     myTree->Branch("RecMuonEta",RecMuonEta,"RecMuonEta[sort_index_for_mu_tree]/F");
967     myTree->Branch("RecMuonPhi",RecMuonPhi,"RecMuonPhi[sort_index_for_mu_tree]/F");
968     myTree->Branch("RecMuonPx",RecMuonPx,"RecMuonPx[sort_index_for_mu_tree]/F");
969     myTree->Branch("RecMuonPy",RecMuonPy,"RecMuonPy[sort_index_for_mu_tree]/F");
970     myTree->Branch("RecMuonPz",RecMuonPz,"RecMuonPz[sort_index_for_mu_tree]/F");
971     myTree->Branch("RecMuonVx",RecMuonVx,"RecMuonVx[sort_index_for_mu_tree]/F");
972     myTree->Branch("RecMuonVy",RecMuonVy,"RecMuonVy[sort_index_for_mu_tree]/F");
973     myTree->Branch("RecMuonVz",RecMuonVz,"RecMuonVz[sort_index_for_mu_tree]/F");
974     myTree->Branch("RecMuonE",RecMuonE,"RecMuonE[sort_index_for_mu_tree]/F");
975     myTree->Branch("RecMuonM",RecMuonM,"RecMuonM[sort_index_for_mu_tree]/F");
976     myTree->Branch("RecMuonGlobalType",RecMuonGlobalType,"RecMuonGlobalType[sort_index_for_mu_tree]/F");
977     myTree->Branch("RecMuonTrackerType",RecMuonTrackerType,"RecMuonTrackerType[sort_index_for_mu_tree]/F");
978     myTree->Branch("RecMuonStandAloneType",RecMuonStandAloneType,"RecMuonStandAloneType[sort_index_for_mu_tree]/F");
979     myTree->Branch("RecMuonIsoSumPt",RecMuonIsoSumPt,"RecMuonIsoSumPt[sort_index_for_mu_tree]/F");
980     myTree->Branch("RecNumberOfUsedStations",RecNumberOfUsedStations,"RecNumberOfUsedStations[sort_index_for_mu_tree]/I");
981     myTree->Branch("RecMuonIsoRelative",RecMuonIsoRelative,"RecMuonIsoRelative[sort_index_for_mu_tree]/F");
982     myTree->Branch("RecMuonIsoCalComb",RecMuonIsoCalComb,"RecMuonIsoCalComb[sort_index_for_mu_tree]/F");
983     myTree->Branch("RecMuonglmuon_dxy",RecMuonglmuon_dxy,"RecMuonglmuon_dxy[sort_index_for_mu_tree]/F");
984     myTree->Branch("RecMuonglmuon_dz",RecMuonglmuon_dz,"RecMuonglmuon_dz[sort_index_for_mu_tree]/F");
985     myTree->Branch("RecMuonglmuon_normalizedChi2",RecMuonglmuon_normalizedChi2,"RecMuonglmuon_normalizedChi2[sort_index_for_mu_tree]/F");
986     myTree->Branch("RecMuonglmuon_trackerHits",RecMuonglmuon_trackerHits,"RecMuonglmuon_trackerHits[sort_index_for_mu_tree]/I");
987     myTree->Branch("RecMuontkmuon_pixelhits",RecMuontkmuon_pixelhits,"RecMuontkmuon_pixelhits[sort_index_for_mu_tree]/I");
988     myTree->Branch("RecMuonglmuon_muonHits",RecMuonglmuon_muonHits,"RecMuonglmuon_muonHits[sort_index_for_mu_tree]/I");
989    
990     myTree->Branch("RecMuonglmuon_charge",RecMuonglmuon_charge,"RecMuonglmuon_charge[sort_index_for_mu_tree]/I");
991     myTree->Branch("hltmatchedmuon",hltmatchedmuon,"hltmatchedmuon[sort_index_for_mu_tree]/I");
992 efe 1.6 myTree->Branch("hltmatched_Dimuon",hltmatched_Dimuon,"hltmatched_Dimuon[sort_index_for_mu_tree]/I");
993     // myTree->Branch("hltmatchedmuon2",hltmatchedmuon2,"hltmatchedmuon2[sort_index_for_mu_tree]/I");
994 efe 1.1
995 efe 1.2 myTree->Branch("recelec_index",&recelec_index,"recelec_index/I");
996     myTree->Branch("RecElec_Pt",RecElec_Pt,"RecElec_Pt[recelec_index]/F");
997     myTree->Branch("RecElec_Px",RecElec_Px,"RecElec_Px[recelec_index]/F");
998     myTree->Branch("RecElec_Py",RecElec_Py,"RecElec_Py[recelec_index]/F");
999     myTree->Branch("RecElec_Pz",RecElec_Pz,"RecElec_Pz[recelec_index]/F");
1000     myTree->Branch("RecElec_eta",RecElec_eta,"RecElec_eta[recelec_index]/F");
1001     myTree->Branch("RecElec_phi",RecElec_phi,"RecElec_phi[recelec_index]/F");
1002     myTree->Branch("RecElec_Charge",RecElec_Charge,"RecElec_Charge[recelec_index]/I");
1003     myTree->Branch("RecElec_GsfTrk_d0",RecElec_GsfTrk_d0,"RecElec_GsfTrk_d0[recelec_index]/F");
1004     myTree->Branch("RecElec_IsEB",RecElec_IsEB,"RecElec_IsEB[recelec_index]/I");
1005     myTree->Branch("RecElec_IsEE",RecElec_IsEE,"RecElec_IsEE[recelec_index]/I");
1006     myTree->Branch("RecElec_dr03TkSumPt",RecElec_dr03TkSumPt,"RecElec_dr03TkSumPt[recelec_index]/F");
1007     myTree->Branch("RecElec_dr03EcalRecHitSumEt",RecElec_dr03EcalRecHitSumEt,"RecElec_dr03EcalRecHitSumEt[recelec_index]/F");
1008     myTree->Branch("RecElec_dr03HcalTowerSumEt",RecElec_dr03HcalTowerSumEt,"RecElec_dr03HcalTowerSumEt[recelec_index]/F");
1009     myTree->Branch("RecElec_scSigmaIEtaIEta",RecElec_scSigmaIEtaIEta,"RecElec_scSigmaIEtaIEta[recelec_index]/F");
1010     myTree->Branch("RecElec_deltaPhiSuperClusterTrackAtVtx",RecElec_deltaPhiSuperClusterTrackAtVtx,"RecElec_deltaPhiSuperClusterTrackAtVtx[recelec_index]/F");
1011     myTree->Branch("RecElec_deltaEtaSuperClusterTrackAtVtx",RecElec_deltaEtaSuperClusterTrackAtVtx,"RecElec_deltaEtaSuperClusterTrackAtVtx[recelec_index]/F");
1012     myTree->Branch("RecElec_hadronicOverEm",RecElec_hadronicOverEm,"RecElec_hadronicOverEm[recelec_index]/F");
1013     myTree->Branch("RecElec_gsfTrack_numberOfLostHits",RecElec_gsfTrack_numberOfLostHits,"RecElec_gsfTrack_numberOfLostHits[recelec_index]/F");
1014    
1015    
1016    
1017 efe 1.1 // myTree->Branch("id_muon",id_muon,"id_muon[sort_index_for_mu_tree]/I");
1018     myTree->Branch("techTrigger",techTrigger, "techTrigger[44]/I");
1019    
1020     //particles
1021    
1022     myTree->Branch("par_index", &par_index, "par_index/I");
1023     myTree->Branch("ParticlePt", ParticlePt, "ParticlePt[par_index]/F");
1024     myTree->Branch("ParticleEta", ParticleEta, "ParticleEta[par_index]/F");
1025     myTree->Branch("ParticlePhi", ParticlePhi, "ParticlePhi[par_index]/F");
1026     myTree->Branch("ParticlePx", ParticlePx, "ParticlePx[par_index]/F");
1027     myTree->Branch("ParticlePy", ParticlePy, "ParticlePy[par_index]/F");
1028     myTree->Branch("ParticlePz", ParticlePz, "ParticlePz[par_index]/F");
1029     myTree->Branch("ParticleE", ParticleE, "ParticleE[par_index]/F");
1030     myTree->Branch("ParticleM", ParticleM, "ParticleM[par_index]/F");
1031     myTree->Branch("ParticleId", ParticleId, "ParticleId[par_index]/I");
1032     myTree->Branch("mom", mom, "mom[par_index]/I");
1033     myTree->Branch("daug", daug, "daug[par_index]/I");
1034     myTree->Branch("ParticleStatus", ParticleStatus, "ParticleStatus[par_index]/I");
1035     myTree->Branch("ParticleMother", ParticleMother, "ParticleMother[par_index][5]/I");
1036     myTree->Branch("ParticleDaughter", ParticleDaughter, "ParticleDaughter[par_index][5]/I");
1037    
1038     //recojets
1039     myTree->Branch("reco_jet",&reco_jet,"reco_jet/I");
1040    
1041     myTree->Branch("RecCorrJetPt",RecCorrJetPt,"RecCorrJetPt[reco_jet]/F");
1042     myTree->Branch("JetEMF",JetEMF,"JetEMF[reco_jet]/F");
1043     myTree->Branch("JetN90",JetN90,"JetN90[reco_jet]/F");
1044     myTree->Branch("JetFHPD",JetFHPD,"JetFHPD[reco_jet]/F");
1045     myTree->Branch("JetFRBX",JetFRBX,"JetFRBX[reco_jet]/F");
1046    
1047     myTree->Branch("RecJetPt",RecJetPt,"RecJetPt[reco_jet]/F");
1048     myTree->Branch("RecJetEta",RecJetEta,"RecJetEta[reco_jet]/F");
1049     myTree->Branch("RecJetPhi",RecJetPhi,"RecJetPhi[reco_jet]/F");
1050     myTree->Branch("RecJetPx",RecJetPx,"RecJetPx[reco_jet]/F");
1051     myTree->Branch("RecJetPy",RecJetPy,"RecJetPy[reco_jet]/F");
1052     myTree->Branch("RecJetPz",RecJetPz,"RecJetPz[reco_jet]/F");
1053     myTree->Branch("RecJetE",RecJetE,"RecJetE[reco_jet]/F");
1054     myTree->Branch("caloMET", &caloMET, "caloMET/F");
1055     myTree->Branch("caloSET", &caloSET, "caloSET/F");
1056     myTree->Branch("pfMET", &pfMET, "pfMET/F");
1057     myTree->Branch("pfSET", &pfSET, "pfSET/F");
1058     myTree->Branch("caloMETX", &caloMETX, "caloMETX/F");
1059     myTree->Branch("pfMETX", &pfMETX, "pfMETX/F");
1060     myTree->Branch("caloMETY", &caloMETY, "caloMETY/F");
1061     myTree->Branch("pfMETY", &pfMETY, "pfMETY/F");
1062     myTree->Branch("muCorrMET", &muCorrMET, "muCorrMET/F");
1063     myTree->Branch("muCorrSET", &muCorrSET, "muCorrSET/F");
1064    
1065     /*
1066     myTree->Branch("jpt_c",&jpt_c,"jpt_c/I");
1067     myTree->Branch("JPTPt",JPTPt,"JPTPt[jpt_c]/F");
1068     myTree->Branch("JPTEta",JPTEta,"JPTEta[jpt_c]/F");
1069     myTree->Branch("JPTPhi",JPTPhi,"JPTPhi[jpt_c]/F");
1070     myTree->Branch("JPTPx",JPTPx,"JPTPx[jpt_c]/F");
1071     myTree->Branch("JPTPy",JPTPy,"JPTPy[jpt_c]/F");
1072     myTree->Branch("JPTPz",JPTPz,"JPTPz[jpt_c]/F");
1073     myTree->Branch("JPTE",JPTE,"JPTE[jpt_c]/F");
1074     */
1075     /*
1076     myTree->Branch("cjpt_c",&cjpt_c,"cjpt_c/I");
1077     myTree->Branch("cJPTPt",cJPTPt,"cJPTPt[cjpt_c]/F");
1078     myTree->Branch("cJPTEta",cJPTEta,"cJPTEta[cjpt_c]/F");
1079     myTree->Branch("cJPTPhi",cJPTPhi,"cJPTPhi[cjpt_c]/F");
1080     myTree->Branch("cJPTPx",cJPTPx,"cJPTPx[cjpt_c]/F");
1081     myTree->Branch("cJPTPy",cJPTPy,"cJPTPy[cjpt_c]/F");
1082     myTree->Branch("cJPTPz",cJPTPz,"cJPTPz[cjpt_c]/F");
1083     myTree->Branch("cJPTE",cJPTE,"cJPTE[cjpt_c]/F");
1084     */
1085    
1086     myTree->Branch("pfNjets",&pfNjets,"pfNjets/I");
1087     myTree->Branch("PFjetEta",PFjetEta,"PFjetEta[pfNjets]/F");
1088     myTree->Branch("PFjetPhi",PFjetPhi,"PFjetPhi[pfNjets]/F");
1089     myTree->Branch("PFjetPt",PFjetPt,"PFjetPt[pfNjets]/F");
1090     myTree->Branch("PFCorrjetPt",PFCorrjetPt,"PFCorrjetPt[pfNjets]/F");
1091     myTree->Branch("PFjetCEMF",PFjetCEMF,"PFjetCEMF[pfNjets]/F");
1092     myTree->Branch("PFjetNEMF",PFjetNEMF,"PFjetNEMF[pfNjets]/F");
1093    
1094     myTree->Branch("pfNtracks",&pfNtracks,"pfNtracks/I");
1095     myTree->Branch("PFjetTrkVZ",PFjetTrkVZ,"PFjetTrkVZ[pfNjets][30]/F");
1096     myTree->Branch("PFjetTrkPT",PFjetTrkPT,"PFjetTrkPT[pfNjets][30]/F");
1097    
1098    
1099     myTree->Branch("nVertices",&nVertices,"nVertices/I");
1100     myTree->Branch("nGoodVertices",&nGoodVertices,"nGoodVertices/I");
1101     myTree->Branch("vtxX",vtxX,"vtxX[nVertices]/F");
1102     myTree->Branch("vtxY",vtxY,"vtxY[nVertices]/F");
1103     myTree->Branch("vtxZ",vtxZ,"vtxZ[nVertices]/F");
1104     myTree->Branch("vtxXerr",vtxXerr,"vtxXerr[nVertices]/F");
1105     myTree->Branch("vtxYerr",vtxYerr,"vtxYerr[nVertices]/F");
1106     myTree->Branch("vtxZerr",vtxZerr,"vtxZerr[nVertices]/F");
1107     myTree->Branch("vtxisValid",vtxisValid,"vtxisValid[nVertices]/I");
1108     myTree->Branch("vtxisFake",vtxisFake,"vtxisFake[nVertices]/I");
1109    
1110     }
1111    
1112     // ------------ method called once each job just after ending the event loop ------------
1113     void
1114     AFB::endJob() {
1115     myTree->Print();//!!!
1116     // f->Write();
1117     // f->Close();
1118     }
1119    
1120    
1121    
1122    
1123     float AFB::DeltaPhi(float phi1, float phi2)
1124     {
1125     float dphi = phi2 - phi1;
1126     if (fabs(dphi) > 3.14) dphi = 6.28 - fabs(dphi);
1127     return dphi;
1128     }
1129    
1130     float AFB::DeltaR(float eta1, float eta2, float phi1, float phi2)
1131     {
1132     float deta = eta2 - eta1;
1133     float dphi = phi2 - phi1;
1134     if (fabs(dphi) > 3.14) dphi = 6.28 - fabs(dphi);
1135     float DELTAR = sqrt(pow(dphi,2)+pow(deta,2))*1.0;
1136     return DELTAR;
1137     }
1138    
1139    
1140     bool AFB::IsMuMatchedToHLTMu ( const reco::Muon & mu, std::vector<reco::Particle> HLTMu , double DR, double DPtRel ) {
1141     size_t dim = HLTMu.size();
1142     size_t nPass=0;
1143    
1144     // filling the denumerator;
1145     double muRecoPt= mu.pt();
1146     //hTrigMuonPtDenS_-> Fill(muRecoPt);
1147    
1148     if (dim==0) return false;
1149     for (size_t k =0; k< dim; k++ ) {
1150     if ( (deltaR(HLTMu[k], mu) < DR) && (fabs(HLTMu[k].pt() - mu.pt())/ HLTMu[k].pt()<DPtRel)){
1151     nPass++ ;
1152     std::cout<< " matching a muon " << std::endl;
1153     std::cout << "muon reco pt : " << muRecoPt<< std::endl;
1154     std::cout << "muon reco eta " << mu.eta() << std::endl;
1155     std::cout << "muon trigger pt "<< HLTMu[k].pt() << std::endl;
1156     std::cout << "muon trigger eta : "<< HLTMu[k].eta() << std::endl;
1157     std::cout <<"deltaR((HLTMu[k], mu)): "<< deltaR(HLTMu[k], mu) << std::endl;
1158     std::cout <<"deltaPtOverPt: "<< fabs(HLTMu[k].pt() - mu.pt())/ HLTMu[k].pt() << std::endl;
1159     }
1160     }
1161    
1162     return (nPass>0);
1163     }
1164    
1165    
1166    
1167     //define this as a plug-in
1168     DEFINE_FWK_MODULE(AFB);