ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc
Revision: 1.18
Committed: Fri May 11 18:43:10 2012 UTC (12 years, 11 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
CVS Tags: HiForest_V02_04
Changes since 1.17: +20 -1 lines
Log Message:
Jet ID variables and calo pf matching

File Contents

# User Rev Content
1 yilmaz 1.1 /*
2     Based on the jet response analyzer
3     Modified by Matt Nguyen, November 2010
4    
5     */
6    
7     #include "CmsHi/JetAnalysis/interface/HiInclusiveJetAnalyzer.h"
8    
9     #include "FWCore/Framework/interface/EventSetup.h"
10     #include "FWCore/Framework/interface/ESHandle.h"
11     #include <Math/DistFunc.h>
12     #include "TMath.h"
13    
14     #include "FWCore/Framework/interface/ESHandle.h"
15    
16     #include "FWCore/MessageLogger/interface/MessageLogger.h"
17     #include "FWCore/Utilities/interface/Exception.h"
18     #include "FWCore/Framework/interface/EventSetup.h"
19    
20     #include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
21     #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
22    
23     #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
24     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
25     #include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
26     #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
27 mnguyen 1.10 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
28 yilmaz 1.1 #include "DataFormats/JetReco/interface/GenJetCollection.h"
29     #include "DataFormats/VertexReco/interface/Vertex.h"
30     #include "DataFormats/Math/interface/deltaPhi.h"
31     #include "DataFormats/Math/interface/deltaR.h"
32    
33     #include "DataFormats/Common/interface/TriggerResults.h"
34     #include "FWCore/Common/interface/TriggerNames.h"
35     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
36     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
37     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
38     #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
39     #include "L1Trigger/GlobalTrigger/interface/L1GlobalTrigger.h"
40    
41     using namespace std;
42     using namespace edm;
43     using namespace reco;
44    
45     HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) {
46    
47     jetTag_ = iConfig.getParameter<InputTag>("jetTag");
48 yilmaz 1.18 matchTag_ = iConfig.getUntrackedParameter<InputTag>("matchTag",jetTag_);
49    
50 mnguyen 1.8 vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));
51 yilmaz 1.17 trackTag_ = iConfig.getParameter<InputTag>("trackTag");
52     useQuality_ = iConfig.getUntrackedParameter<bool>("useQuality",0);
53     string trackQuality_ = iConfig.getUntrackedParameter<string>("trackQuality","highPurity");
54 yilmaz 1.1
55     isMC_ = iConfig.getUntrackedParameter<bool>("isMC",false);
56 yilmaz 1.14 doTrigger_ = iConfig.getUntrackedParameter<bool>("doTrigger",false);
57 yilmaz 1.17
58     rParam = iConfig.getParameter<double>("rParam");
59 mnguyen 1.8
60 yilmaz 1.1 if(isMC_){
61     genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
62     eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
63     }
64     verbose_ = iConfig.getUntrackedParameter<bool>("verbose",false);
65    
66     useCentrality_ = iConfig.getUntrackedParameter<bool>("useCentrality",false);
67 yjlee 1.4 useVtx_ = iConfig.getUntrackedParameter<bool>("useVtx",false);
68 yilmaz 1.1 useJEC_ = iConfig.getUntrackedParameter<bool>("useJEC",true);
69 yilmaz 1.2 usePat_ = iConfig.getUntrackedParameter<bool>("usePAT",true);
70 yilmaz 1.1
71 mnguyen 1.8 doLifeTimeTagging_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTagging",false);
72 mnguyen 1.10 doLifeTimeTaggingExtras_ = iConfig.getUntrackedParameter<bool>("doLifeTimeTaggingExtras",true);
73 mnguyen 1.11 saveBfragments_ = iConfig.getUntrackedParameter<bool>("saveBfragments",false);
74    
75 mnguyen 1.8 pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
76    
77 yilmaz 1.14 if(doTrigger_){
78 yilmaz 1.1 L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
79     hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
80    
81    
82     if (iConfig.exists("hltTrgNames"))
83     hltTrgNames_ = iConfig.getUntrackedParameter<vector<string> >("hltTrgNames");
84    
85     if (iConfig.exists("hltProcNames"))
86     hltProcNames_ = iConfig.getUntrackedParameter<vector<string> >("hltProcNames");
87     else {
88     hltProcNames_.push_back("FU");
89     hltProcNames_.push_back("HLT");
90     }
91     }
92    
93     cout<<" jet collection : "<<jetTag_<<endl;
94 yilmaz 1.12 doSubEvent_ = 0;
95    
96 yilmaz 1.9 if(isMC_){
97     cout<<" genjet collection : "<<genjetTag_<<endl;
98     genPtMin_ = iConfig.getUntrackedParameter<double>("genPtMin",0);
99 yilmaz 1.12 doSubEvent_ = iConfig.getUntrackedParameter<bool>("doSubEvent",1);
100 yilmaz 1.9 }
101 yilmaz 1.1
102    
103     }
104    
105    
106    
107     HiInclusiveJetAnalyzer::~HiInclusiveJetAnalyzer() { }
108    
109    
110    
111     void
112     HiInclusiveJetAnalyzer::beginRun(const edm::Run& run,
113     const edm::EventSetup & es) {}
114    
115     void
116     HiInclusiveJetAnalyzer::beginJob() {
117    
118     centrality_ = 0;
119    
120     //string jetTagName = jetTag_.label()+"_tree";
121     string jetTagTitle = jetTag_.label()+" Jet Analysis Tree";
122     t = fs1->make<TTree>("t",jetTagTitle.c_str());
123    
124 yjlee 1.6 // TTree* t= new TTree("t","Jet Response Analyzer");
125 yjlee 1.4 //t->Branch("run",&jets_.run,"run/I");
126 yilmaz 1.1 t->Branch("evt",&jets_.evt,"evt/I");
127 yjlee 1.4 //t->Branch("lumi",&jets_.lumi,"lumi/I");
128 yilmaz 1.1 t->Branch("b",&jets_.b,"b/F");
129 yjlee 1.4 if (useVtx_) {
130     t->Branch("vx",&jets_.vx,"vx/F");
131     t->Branch("vy",&jets_.vy,"vy/F");
132     t->Branch("vz",&jets_.vz,"vz/F");
133     }
134     if (useCentrality_) {
135     t->Branch("hf",&jets_.hf,"hf/F");
136     t->Branch("bin",&jets_.bin,"bin/I");
137     }
138 yjlee 1.7
139     t->Branch("nref",&jets_.nref,"nref/I");
140 yilmaz 1.1 t->Branch("rawpt",jets_.rawpt,"rawpt[nref]/F");
141     t->Branch("jtpt",jets_.jtpt,"jtpt[nref]/F");
142     t->Branch("jteta",jets_.jteta,"jteta[nref]/F");
143     t->Branch("jty",jets_.jty,"jty[nref]/F");
144     t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
145 yilmaz 1.2 t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
146 yilmaz 1.1
147 yilmaz 1.17 // jet ID information, jet composition
148     t->Branch("discr_fr01", jets_.discr_fr01,"discr_fr01[nref]/F");
149    
150     t->Branch("trackMax", jets_.trackMax,"trackMax[nref]/F");
151     t->Branch("trackSum", jets_.trackSum,"trackSum[nref]/F");
152     t->Branch("trackN", jets_.trackN,"trackN[nref]/I");
153    
154     t->Branch("chargedMax", jets_.chargedMax,"chargedMax[nref]/F");
155     t->Branch("chargedSum", jets_.chargedSum,"chargedSum[nref]/F");
156     t->Branch("chargedN", jets_.chargedN,"chargedN[nref]/I");
157    
158     t->Branch("photonMax", jets_.photonMax,"photonMax[nref]/F");
159     t->Branch("photonSum", jets_.photonSum,"photonSum[nref]/F");
160     t->Branch("photonN", jets_.photonN,"photonN[nref]/I");
161    
162     t->Branch("neutralMax", jets_.neutralMax,"neutralMax[nref]/F");
163     t->Branch("neutralSum", jets_.neutralSum,"neutralSum[nref]/F");
164     t->Branch("neutralN", jets_.neutralN,"neutralN[nref]/I");
165    
166     t->Branch("eMax", jets_.eMax,"eMax[nref]/F");
167     t->Branch("eSum", jets_.eSum,"eSum[nref]/F");
168     t->Branch("eN", jets_.eN,"eN[nref]/I");
169    
170     t->Branch("muMax", jets_.muMax,"muMax[nref]/F");
171     t->Branch("muSum", jets_.muSum,"muSum[nref]/F");
172     t->Branch("muN", jets_.muN,"muN[nref]/I");
173    
174 yilmaz 1.18 t->Branch("matchedPt", jets_.matchedPt,"matchedPt[nref]/F");
175     t->Branch("matchedR", jets_.matchedR,"matchedR[nref]/F");
176 yilmaz 1.17
177 mnguyen 1.8 // b-jet discriminators
178     if (doLifeTimeTagging_) {
179    
180     t->Branch("discr_csvMva",jets_.discr_csvMva,"discr_csvMva[nref]/F");
181     t->Branch("discr_csvSimple",jets_.discr_csvSimple,"discr_csvSimple[nref]/F");
182     t->Branch("discr_muByIp3",jets_.discr_muByIp3,"discr_muByIp3[nref]/F");
183     t->Branch("discr_muByPt",jets_.discr_muByPt,"discr_muByPt[nref]/F");
184     t->Branch("discr_prob",jets_.discr_prob,"discr_prob[nref]/F");
185     t->Branch("discr_probb",jets_.discr_probb,"discr_probb[nref]/F");
186     t->Branch("discr_tcHighEff",jets_.discr_tcHighEff,"discr_tcHighEff[nref]/F");
187     t->Branch("discr_tcHighPur",jets_.discr_tcHighPur,"discr_tcHighPur[nref]/F");
188    
189 mnguyen 1.10 t->Branch("nsvtx", jets_.nsvtx, "nsvtx[nref]/I");
190     t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/I");
191 mnguyen 1.8 t->Branch("svtxdl", jets_.svtxdl, "svtxdl[nref]/F");
192     t->Branch("svtxdls", jets_.svtxdls, "svtxdls[nref]/F");
193     t->Branch("svtxm", jets_.svtxm, "svtxm[nref]/F");
194     t->Branch("svtxpt", jets_.svtxpt, "svtxpt[nref]/F");
195    
196 mnguyen 1.10 t->Branch("nIPtrk",jets_.nIPtrk,"nIPtrk[nref]/I");
197     t->Branch("nselIPtrk",jets_.nselIPtrk,"nselIPtrk[nref]/I");
198 mnguyen 1.11
199 mnguyen 1.10 if (doLifeTimeTaggingExtras_) {
200 mnguyen 1.11 t->Branch("nIP",&jets_.nIP,"nIP/I");
201     t->Branch("ipJetIndex",jets_.ipJetIndex,"ipJetIndex[nIP]/I");
202     t->Branch("ipPt",jets_.ipPt,"ipPt[nIP]/F");
203     t->Branch("ipProb0",jets_.ipProb0,"ipProb0[nIP]/F");
204     t->Branch("ipProb1",jets_.ipProb1,"ipProb1[nIP]/F");
205     t->Branch("ip2d",jets_.ip2d,"ip2d[nIP]/F");
206     t->Branch("ip2dSig",jets_.ip2dSig,"ip2dSig[nIP]/F");
207     t->Branch("ip3d",jets_.ip3d,"ip3d[nIP]/F");
208     t->Branch("ip3dSig",jets_.ip3dSig,"ip3dSig[nIP]/F");
209     t->Branch("ipDist2Jet",jets_.ipDist2Jet,"ipDist2Jet[nIP]/F");
210     t->Branch("ipDist2JetSig",jets_.ipDist2JetSig,"ipDist2JetSig[nIP]/F");
211     t->Branch("ipClosest2Jet",jets_.ipClosest2Jet,"ipClosest2Jet[nIP]/F");
212    
213 mnguyen 1.10 }
214 mnguyen 1.8
215     t->Branch("mue", jets_.mue, "mue[nref]/F");
216     t->Branch("mupt", jets_.mupt, "mupt[nref]/F");
217     t->Branch("mueta", jets_.mueta, "mueta[nref]/F");
218     t->Branch("muphi", jets_.muphi, "muphi[nref]/F");
219     t->Branch("mudr", jets_.mudr, "mudr[nref]/F");
220     t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
221     t->Branch("muchg", jets_.muchg, "muchg[nref]/I");
222     }
223 ylai 1.16
224 mnguyen 1.8
225 yilmaz 1.1 if(isMC_){
226 mnguyen 1.10 t->Branch("beamId1",&jets_.beamId1,"beamId1/I");
227     t->Branch("beamId2",&jets_.beamId2,"beamId2/I");
228    
229 yilmaz 1.1 t->Branch("pthat",&jets_.pthat,"pthat/F");
230    
231     // Only matched gen jets
232     t->Branch("refpt",jets_.refpt,"refpt[nref]/F");
233     t->Branch("refeta",jets_.refeta,"refeta[nref]/F");
234     t->Branch("refy",jets_.refy,"refy[nref]/F");
235     t->Branch("refphi",jets_.refphi,"refphi[nref]/F");
236     t->Branch("refdphijt",jets_.refdphijt,"refdphijt[nref]/F");
237     t->Branch("refdrjt",jets_.refdrjt,"refdrjt[nref]/F");
238     // matched parton
239     t->Branch("refparton_pt",jets_.refparton_pt,"refparton_pt[nref]/F");
240 mnguyen 1.8 t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
241     t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
242 yilmaz 1.1
243     // For all gen jets, matched or unmatched
244     t->Branch("ngen",&jets_.ngen,"ngen/I");
245     t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
246     t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
247     t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
248     t->Branch("geny",jets_.geny,"geny[ngen]/F");
249     t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
250     t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
251     t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
252 yilmaz 1.12
253     if(doSubEvent_){
254     t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
255     }
256 mnguyen 1.11
257     if(saveBfragments_ ) {
258     t->Branch("bMult",&jets_.bMult,"bMult/I");
259     t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
260     t->Branch("bStatus",jets_.bStatus,"bStatus[bMult]/I");
261     t->Branch("bVx",jets_.bVx,"bVx[bMult]/F");
262     t->Branch("bVy",jets_.bVy,"bVy[bMult]/F");
263     t->Branch("bVz",jets_.bVz,"bVz[bMult]/F");
264     t->Branch("bPt",jets_.bPt,"bPt[bMult]/F");
265     t->Branch("bEta",jets_.bEta,"bEta[bMult]/F");
266     t->Branch("bPhi",jets_.bPhi,"bPhi[bMult]/F");
267     t->Branch("bPdg",jets_.bPdg,"bPdg[bMult]/I");
268     t->Branch("bChg",jets_.bChg,"bChg[bMult]/I");
269     }
270 yilmaz 1.13
271 yilmaz 1.1 }
272 yjlee 1.7 /*
273 yilmaz 1.1 if(!isMC_){
274     t->Branch("nL1TBit",&jets_.nL1TBit,"nL1TBit/I");
275     t->Branch("l1TBit",jets_.l1TBit,"l1TBit[nL1TBit]/O");
276    
277     t->Branch("nL1ABit",&jets_.nL1ABit,"nL1ABit/I");
278     t->Branch("l1ABit",jets_.l1ABit,"l1ABit[nL1ABit]/O");
279    
280     t->Branch("nHLTBit",&jets_.nHLTBit,"nHLTBit/I");
281     t->Branch("hltBit",jets_.hltBit,"hltBit[nHLTBit]/O");
282    
283     }
284 yjlee 1.7 */
285 yilmaz 1.1 TH1D::SetDefaultSumw2();
286    
287    
288     }
289    
290    
291     void
292     HiInclusiveJetAnalyzer::analyze(const Event& iEvent,
293     const EventSetup& iSetup) {
294    
295     int event = iEvent.id().event();
296     int run = iEvent.id().run();
297     int lumi = iEvent.id().luminosityBlock();
298    
299     jets_.run = run;
300     jets_.evt = event;
301     jets_.lumi = lumi;
302    
303     LogDebug("HiInclusiveJetAnalyzer")<<"START event: "<<event<<" in run "<<run<<endl;
304    
305     int bin = -1;
306     double hf = 0.;
307     double b = 999.;
308    
309     if(useCentrality_){
310     if(!centrality_) centrality_ = new CentralityProvider(iSetup);
311     centrality_->newEvent(iEvent,iSetup); // make sure you do this first in every event
312     //double c = centrality_->centralityValue();
313     const reco::Centrality *cent = centrality_->raw();
314    
315     hf = cent->EtHFhitSum();
316    
317     bin = centrality_->getBin();
318     b = centrality_->bMean();
319     }
320    
321    
322    
323    
324    
325     // loop the events
326    
327     jets_.bin = bin;
328     jets_.hf = hf;
329    
330    
331     if (useVtx_) {
332     edm::Handle<vector<reco::Vertex> >vertex;
333     iEvent.getByLabel(vtxTag_, vertex);
334    
335     if(vertex->size()>0) {
336     jets_.vx=vertex->begin()->x();
337     jets_.vy=vertex->begin()->y();
338     jets_.vz=vertex->begin()->z();
339     }
340     }
341    
342 yilmaz 1.2 edm::Handle<pat::JetCollection> patjets;
343     if(usePat_)iEvent.getByLabel(jetTag_, patjets);
344 yilmaz 1.18
345     edm::Handle<pat::JetCollection> matchedjets;
346     iEvent.getByLabel(matchTag_, matchedjets);
347 yilmaz 1.2
348     edm::Handle<reco::JetView> jets;
349 yilmaz 1.1 iEvent.getByLabel(jetTag_, jets);
350    
351 mnguyen 1.8 edm::Handle<reco::PFCandidateCollection> pfCandidates;
352 yilmaz 1.17 iEvent.getByLabel(pfCandidateLabel_,pfCandidates);
353    
354     edm::Handle<reco::TrackCollection> tracks;
355     iEvent.getByLabel(trackTag_,tracks);
356    
357 yilmaz 1.1 // FILL JRA TREE
358    
359     jets_.b = b;
360     jets_.nref = 0;
361    
362 yilmaz 1.15 if(doTrigger_){
363 yilmaz 1.1 fillL1Bits(iEvent);
364     fillHLTBits(iEvent);
365     }
366    
367     for(unsigned int j = 0 ; j < jets->size(); ++j){
368 yilmaz 1.2 const reco::Jet& jet = (*jets)[j];
369 yilmaz 1.1
370 mnguyen 1.11
371 yilmaz 1.2 if (useJEC_ && usePat_){
372     jets_.rawpt[jets_.nref]=(*patjets)[j].correctedJet("Uncorrected").pt();
373     }
374 mnguyen 1.8
375     if(doLifeTimeTagging_){
376    
377     if(jetTag_.label()=="icPu5patJets"){
378     jets_.discr_csvMva[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexMVABJetTags");
379     jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5CombinedSecondaryVertexBJetTags");
380     jets_.discr_muByIp3[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5SoftMuonByIP3dBJetTags");
381     jets_.discr_muByPt[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5SoftMuonByPtBJetTags");
382     jets_.discr_prob[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5JetBProbabilityBJetTags");
383     jets_.discr_probb[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5JetProbabilityBJetTags");
384     jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighEffBJetTags");
385     jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("icPu5TrackCountingHighPurBJetTags");
386     }
387     else if(jetTag_.label()=="akPu3PFpatJets"){
388     jets_.discr_csvMva[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexMVABJetTags");
389     jets_.discr_csvSimple[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFCombinedSecondaryVertexBJetTags");
390     jets_.discr_muByIp3[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByIP3dBJetTags");
391     jets_.discr_muByPt[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFSoftMuonByPtBJetTags");
392     jets_.discr_prob[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFJetBProbabilityBJetTags");
393     jets_.discr_probb[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFJetProbabilityBJetTags");
394     jets_.discr_tcHighEff[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighEffBJetTags");
395     jets_.discr_tcHighPur[jets_.nref] = (*patjets)[j].bDiscriminator("akPu3PFTrackCountingHighPurBJetTags");
396     }
397     else{
398 mnguyen 1.10 cout<<" b-tagging variables not filled for this collection, turn of doLifeTimeTagging "<<endl;
399 mnguyen 1.8 }
400 mnguyen 1.10
401     const reco::SecondaryVertexTagInfo& tagInfoSV=*(*patjets)[j].tagInfoSecondaryVertex();
402    
403     jets_.nsvtx[jets_.nref] = tagInfoSV.nVertices();
404 mnguyen 1.8
405     if (tagInfoSV.nVertices()>0) {
406     jets_.svtxntrk[jets_.nref] = tagInfoSV.nVertexTracks(0);
407 mnguyen 1.10 // this is the 3d flight distance, for 2-D use (0,true)
408     Measurement1D m1D = tagInfoSV.flightDistance(0);
409 mnguyen 1.8 jets_.svtxdl[jets_.nref] = m1D.value();
410     jets_.svtxdls[jets_.nref] = m1D.significance();
411    
412 mnguyen 1.10 const reco::Vertex& svtx = tagInfoSV.secondaryVertex(0);
413 mnguyen 1.11 //cout<<" SV: vx: "<<svtx.x()<<" vy "<<svtx.y()<<" vz "<<svtx.z()<<endl;
414     //cout<<" PV: vx: "<<jet.vx()<<" vy "<<jet.vy()<<" vz "<<jet.vz()<<endl;
415 mnguyen 1.8 jets_.svtxm[jets_.nref] = svtx.p4().mass();
416 mnguyen 1.10 jets_.svtxpt[jets_.nref] = svtx.p4().pt();
417 mnguyen 1.11 //cout<<" chi2 "<<svtx.chi2()<<" ndof "<<svtx.ndof()<<endl;
418 mnguyen 1.8 }
419 mnguyen 1.10
420 mnguyen 1.8 const reco::TrackIPTagInfo& tagInfoIP=*(*patjets)[j].tagInfoTrackIP();
421 mnguyen 1.10
422     jets_.nIPtrk[jets_.nref] = tagInfoIP.tracks().size();
423     jets_.nselIPtrk[jets_.nref] = tagInfoIP.selectedTracks().size();
424 mnguyen 1.8
425 mnguyen 1.10 if (doLifeTimeTaggingExtras_) {
426 mnguyen 1.8
427 mnguyen 1.10 TrackRefVector selTracks=tagInfoIP.selectedTracks();
428    
429     GlobalPoint pv(tagInfoIP.primaryVertex()->position().x(),tagInfoIP.primaryVertex()->position().y(),tagInfoIP.primaryVertex()->position().z());
430    
431     for(int it=0;it<jets_.nselIPtrk[jets_.nref] ;it++)
432     {
433 mnguyen 1.11 jets_.ipJetIndex[jets_.nIP + it]= jets_.nref;
434 mnguyen 1.10 TrackIPTagInfo::TrackIPData data = tagInfoIP.impactParameterData()[it];
435 mnguyen 1.11 jets_.ipPt[jets_.nIP + it] = selTracks[it]->pt();
436     jets_.ipProb0[jets_.nIP + it] = tagInfoIP.probabilities(0)[it];
437     jets_.ipProb1[jets_.nIP + it] = tagInfoIP.probabilities(1)[it];
438     jets_.ip2d[jets_.nIP + it] = data.ip2d.value();
439     jets_.ip2dSig[jets_.nIP + it] = data.ip2d.significance();
440     jets_.ip3d[jets_.nIP + it] = data.ip3d.value();
441     jets_.ip3dSig[jets_.nIP + it] = data.ip3d.significance();
442     jets_.ipDist2Jet[jets_.nIP + it] = data.distanceToJetAxis.value();
443     jets_.ipDist2JetSig[jets_.nIP + it] = data.distanceToJetAxis.significance();
444     jets_.ipClosest2Jet[jets_.nIP + it] = (data.closestToJetAxis - pv).mag(); //decay length
445 mnguyen 1.10 }
446 mnguyen 1.11
447     jets_.nIP += jets_.nselIPtrk[jets_.nref];
448    
449 mnguyen 1.10 }
450    
451 mnguyen 1.8 const reco::PFCandidateCollection *pfCandidateColl = &(*pfCandidates);
452     int pfMuonIndex = getPFJetMuon(jet, pfCandidateColl);
453     if(pfMuonIndex >=0){
454     const reco::PFCandidate muon = pfCandidateColl->at(pfMuonIndex);
455     jets_.mupt[jets_.nref] = muon.pt();
456     jets_.mueta[jets_.nref] = muon.eta();
457     jets_.muphi[jets_.nref] = muon.phi();
458     jets_.mue[jets_.nref] = muon.energy();
459     jets_.mudr[jets_.nref] = reco::deltaR(jet,muon);
460     jets_.muptrel[jets_.nref] = getPtRel(muon, jet);
461     jets_.muchg[jets_.nref] = muon.charge();
462 yilmaz 1.17 }else{
463 ylai 1.16 jets_.mupt[jets_.nref] = 0.0;
464     jets_.mueta[jets_.nref] = 0.0;
465     jets_.muphi[jets_.nref] = 0.0;
466     jets_.mue[jets_.nref] = 0.0;
467     jets_.mudr[jets_.nref] = 9.9;
468     jets_.muptrel[jets_.nref] = 0.0;
469     jets_.muchg[jets_.nref] = 0;
470     }
471    
472 mnguyen 1.8 }
473    
474 yilmaz 1.17
475    
476     // Jet ID variables
477    
478     jets_.muMax[jets_.nref] = 0;
479     jets_.muSum[jets_.nref] = 0;
480     jets_.muN[jets_.nref] = 0;
481    
482     jets_.eMax[jets_.nref] = 0;
483     jets_.eSum[jets_.nref] = 0;
484     jets_.eN[jets_.nref] = 0;
485    
486     jets_.neutralMax[jets_.nref] = 0;
487     jets_.neutralSum[jets_.nref] = 0;
488     jets_.neutralN[jets_.nref] = 0;
489    
490     jets_.photonMax[jets_.nref] = 0;
491     jets_.photonSum[jets_.nref] = 0;
492     jets_.photonN[jets_.nref] = 0;
493    
494     jets_.chargedMax[jets_.nref] = 0;
495     jets_.chargedSum[jets_.nref] = 0;
496     jets_.chargedN[jets_.nref] = 0;
497    
498     jets_.trackMax[jets_.nref] = 0;
499     jets_.trackSum[jets_.nref] = 0;
500     jets_.trackN[jets_.nref] = 0;
501    
502    
503     for(unsigned int icand = 0; icand < tracks->size(); ++icand){
504     const reco::Track& track = (*tracks)[icand];
505 yilmaz 1.18 if(useQuality_ && !(track.quality(reco::TrackBase::qualityByName(trackQuality_)))) continue;
506    
507 yilmaz 1.17 double dr = deltaR(jet,track);
508     if(dr < rParam){
509     double ptcand = track.pt();
510     jets_.trackSum[jets_.nref] += ptcand;
511     jets_.trackN[jets_.nref] += 1;
512     if(ptcand > jets_.trackMax[jets_.nref]) jets_.trackMax[jets_.nref] = ptcand;
513    
514     }
515     }
516    
517     for(unsigned int icand = 0; icand < pfCandidates->size(); ++icand){
518     const reco::PFCandidate& track = (*pfCandidates)[icand];
519     double dr = deltaR(jet,track);
520     if(dr < rParam){
521     double ptcand = track.pt();
522     int pfid = track.particleId();
523    
524     switch(pfid){
525    
526     case 1:
527     jets_.chargedSum[jets_.nref] += ptcand;
528     jets_.chargedN[jets_.nref] += 1;
529     if(ptcand > jets_.chargedMax[jets_.nref]) jets_.chargedMax[jets_.nref] = ptcand;
530     break;
531    
532     case 2:
533     jets_.eSum[jets_.nref] += ptcand;
534     jets_.eN[jets_.nref] += 1;
535     if(ptcand > jets_.eMax[jets_.nref]) jets_.eMax[jets_.nref] = ptcand;
536     break;
537    
538     case 3:
539     jets_.muSum[jets_.nref] += ptcand;
540     jets_.muN[jets_.nref] += 1;
541     if(ptcand > jets_.muMax[jets_.nref]) jets_.muMax[jets_.nref] = ptcand;
542     break;
543    
544     case 4:
545     jets_.photonSum[jets_.nref] += ptcand;
546     jets_.photonN[jets_.nref] += 1;
547     if(ptcand > jets_.photonMax[jets_.nref]) jets_.photonMax[jets_.nref] = ptcand;
548     break;
549    
550     case 5:
551     jets_.neutralSum[jets_.nref] += ptcand;
552     jets_.neutralN[jets_.nref] += 1;
553     if(ptcand > jets_.neutralMax[jets_.nref]) jets_.neutralMax[jets_.nref] = ptcand;
554     break;
555    
556     }
557     }
558     }
559    
560 yilmaz 1.18 double drMin = 100;
561     for(unsigned int j = 0 ; j < matchedjets->size(); ++j){
562     const reco::Jet& mjet = (*matchedjets)[j];
563    
564     double dr = deltaR(jet,mjet);
565     if(dr < drMin){
566     jets_.matchedPt[jets_.nref] = mjet.pt();
567     jets_.matchedR[jets_.nref] = dr;
568    
569     }
570     }
571 yilmaz 1.17
572    
573    
574    
575     // if(etrk.quality(reco::TrackBase::qualityByName(qualityString_))) pev_.trkQual[pev_.nTrk]=1;
576    
577 ylai 1.16 /////////////////////////////////////////////////////////////////
578     // Jet core pt^2 discriminant for fake jets
579     // Edited by Yue Shi Lai <ylai@mit.edu>
580    
581     // Initial value is 0
582     jets_.discr_fr01[jets_.nref] = 0;
583     // Start with no directional adaption, i.e. the fake rejection
584     // axis is the jet axis
585     float pseudorapidity_adapt = jets_.jteta[jets_.nref];
586     float azimuth_adapt = jets_.jtphi[jets_.nref];
587    
588     // Unadapted discriminant with adaption search
589     for (size_t iteration = 0; iteration < 2; iteration++) {
590     float pseudorapidity_adapt_new = pseudorapidity_adapt;
591     float azimuth_adapt_new = azimuth_adapt;
592     float max_weighted_perp = 0;
593     float perp_square_sum = 0;
594    
595     for (size_t index_pf_candidate = 0;
596     index_pf_candidate < pfCandidates->size();
597     index_pf_candidate++) {
598     const reco::PFCandidate &p =
599     (*pfCandidates)[index_pf_candidate];
600    
601     switch (p.particleId()) {
602     case 1: // Charged hadron
603     case 3: // Muon
604     case 4: // Photon
605     {
606     const float dpseudorapidity =
607     p.eta() - pseudorapidity_adapt;
608     const float dazimuth =
609     reco::deltaPhi(p.phi(), azimuth_adapt);
610     // The Gaussian scale factor is 0.5 / (0.1 * 0.1)
611     // = 50
612     const float angular_weight =
613     exp(-50.0F * (dpseudorapidity * dpseudorapidity +
614     dazimuth * dazimuth));
615     const float weighted_perp =
616     angular_weight * p.pt() * p.pt();
617     const float weighted_perp_square =
618     weighted_perp * p.pt();
619    
620     perp_square_sum += weighted_perp_square;
621     if (weighted_perp >= max_weighted_perp) {
622     pseudorapidity_adapt_new = p.eta();
623     azimuth_adapt_new = p.phi();
624     max_weighted_perp = weighted_perp;
625     }
626     }
627     }
628     }
629     // Update the fake rejection value
630     jets_.discr_fr01[jets_.nref] = std::max(
631     jets_.discr_fr01[jets_.nref], perp_square_sum);
632     // Update the directional adaption
633     pseudorapidity_adapt = pseudorapidity_adapt_new;
634     azimuth_adapt = azimuth_adapt_new;
635     }
636    
637    
638 yilmaz 1.1 jets_.jtpt[jets_.nref] = jet.pt();
639     jets_.jteta[jets_.nref] = jet.eta();
640     jets_.jtphi[jets_.nref] = jet.phi();
641     jets_.jty[jets_.nref] = jet.eta();
642 yilmaz 1.2 jets_.jtpu[jets_.nref] = jet.pileup();
643 yilmaz 1.1
644 yilmaz 1.2 if(isMC_ && usePat_){
645 mnguyen 1.11
646     const reco::GenJet * genjet = (*patjets)[j].genJet();
647 mnguyen 1.8
648 yilmaz 1.2 if(genjet){
649     jets_.refpt[jets_.nref] = genjet->pt();
650     jets_.refeta[jets_.nref] = genjet->eta();
651     jets_.refphi[jets_.nref] = genjet->phi();
652     jets_.refy[jets_.nref] = genjet->eta();
653     jets_.refdphijt[jets_.nref] = reco::deltaPhi(jet.phi(), genjet->phi());
654     jets_.refdrjt[jets_.nref] = reco::deltaR(jet.eta(),jet.phi(),genjet->eta(),genjet->phi());
655     }else{
656     jets_.refpt[jets_.nref] = -999.;
657     jets_.refeta[jets_.nref] = -999.;
658     jets_.refphi[jets_.nref] = -999.;
659     jets_.refy[jets_.nref] = -999.;
660     jets_.refdphijt[jets_.nref] = -999.;
661     jets_.refdrjt[jets_.nref] = -999.;
662     }
663 yilmaz 1.1
664 mnguyen 1.11 jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour();
665    
666     // matched partons
667     const reco::GenParticle & parton = *(*patjets)[j].genParton();
668    
669     if((*patjets)[j].genParton()){
670     jets_.refparton_pt[jets_.nref] = parton.pt();
671     jets_.refparton_flavor[jets_.nref] = parton.pdgId();
672    
673     if(saveBfragments_ && abs(jets_.refparton_flavorForB[jets_.nref])==5){
674    
675     usedStringPts.clear();
676    
677     // uncomment this if you want to know the ugly truth about parton matching -matt
678     //if(jet.pt() > 50 &&abs(parton.pdgId())!=5 && parton.pdgId()!=21)
679     // cout<<" Identified as a b, but doesn't match b or gluon, id = "<<parton.pdgId()<<endl;
680    
681     jets_.bJetIndex[jets_.bMult] = jets_.nref;
682     jets_.bStatus[jets_.bMult] = parton.status();
683     jets_.bVx[jets_.bMult] = parton.vx();
684     jets_.bVy[jets_.bMult] = parton.vy();
685     jets_.bVz[jets_.bMult] = parton.vz();
686     jets_.bPt[jets_.bMult] = parton.pt();
687     jets_.bEta[jets_.bMult] = parton.eta();
688     jets_.bPhi[jets_.bMult] = parton.phi();
689     jets_.bPdg[jets_.bMult] = parton.pdgId();
690     jets_.bChg[jets_.bMult] = parton.charge();
691     jets_.bMult++;
692     saveDaughters(parton);
693     }
694    
695    
696 mnguyen 1.8 } else {
697 mnguyen 1.11 jets_.refparton_pt[jets_.nref] = -999;
698     jets_.refparton_flavor[jets_.nref] = -999;
699 mnguyen 1.8 }
700 mnguyen 1.11
701    
702 yilmaz 1.2 }
703 yilmaz 1.1
704     jets_.nref++;
705    
706    
707     }
708    
709    
710     if(isMC_){
711    
712 mnguyen 1.10 edm::Handle<HepMCProduct> hepMCProduct;
713     iEvent.getByLabel(eventInfoTag_,hepMCProduct);
714     const HepMC::GenEvent* MCEvt = hepMCProduct->GetEvent();
715     std::pair<HepMC::GenParticle*,HepMC::GenParticle*> beamParticles = MCEvt->beam_particles();
716     jets_.beamId1 = beamParticles.first->pdg_id();
717     jets_.beamId2 = beamParticles.second->pdg_id();
718    
719 yilmaz 1.1 edm::Handle<GenEventInfoProduct> hEventInfo;
720     iEvent.getByLabel(eventInfoTag_,hEventInfo);
721     //jets_.pthat = hEventInfo->binningValues()[0];
722    
723     // binning values and qscale appear to be equivalent, but binning values not always present
724     jets_.pthat = hEventInfo->qScale();
725    
726     edm::Handle<vector<reco::GenJet> >genjets;
727     iEvent.getByLabel(genjetTag_, genjets);
728    
729     jets_.ngen = 0;
730     for(unsigned int igen = 0 ; igen < genjets->size(); ++igen){
731     const reco::GenJet & genjet = (*genjets)[igen];
732    
733     float genjet_pt = genjet.pt();
734    
735     // threshold to reduce size of output in minbias PbPb
736 yilmaz 1.9 if(genjet_pt>genPtMin_){
737 yilmaz 1.1
738     jets_.genpt[jets_.ngen] = genjet_pt;
739     jets_.geneta[jets_.ngen] = genjet.eta();
740     jets_.genphi[jets_.ngen] = genjet.phi();
741     jets_.geny[jets_.ngen] = genjet.eta();
742    
743 yilmaz 1.12 if(doSubEvent_){
744     const GenParticle* gencon = genjet.getGenConstituent(0);
745     jets_.gensubid[jets_.ngen] = gencon->collisionId();
746     }
747 yilmaz 1.1
748     // find matching patJet if there is one
749    
750     jets_.gendrjt[jets_.ngen] = -1.0;
751     jets_.genmatchindex[jets_.ngen] = -1;
752    
753     for(unsigned int ijet = 0 ; ijet < jets->size(); ++ijet){
754     const pat::Jet& jet = (*jets)[ijet];
755 mnguyen 1.8 const reco::GenJet* matchedGenJet = (*patjets)[ijet].genJet();
756     if(matchedGenJet){
757     // poor man's matching, someone fix please
758     if(fabs(genjet.pt()-matchedGenJet->pt())<0.0001 &&
759     fabs(genjet.eta()-matchedGenJet->eta())<0.0001 &&
760     (fabs(genjet.phi()-matchedGenJet->phi())<0.0001 || fabs(genjet.phi()-matchedGenJet->phi() - 2.0*TMath::Pi()) < 0.0001 )){
761 yilmaz 1.1
762     jets_.genmatchindex[jets_.ngen] = (int)ijet;
763     jets_.gendphijt[jets_.ngen] = reco::deltaPhi(jet.phi(),genjet.phi());
764     jets_.gendrjt[jets_.ngen] = reco::deltaR(jet,genjet);
765    
766     break;
767     }
768     }
769     }
770     jets_.ngen++;
771     }
772    
773     }
774    
775     }
776     t->Fill();
777 mnguyen 1.10 memset(&jets_,0,sizeof jets_);
778 yilmaz 1.1 }
779    
780    
781    
782    
783     //--------------------------------------------------------------------------------------------------
784     void HiInclusiveJetAnalyzer::fillL1Bits(const edm::Event &iEvent)
785     {
786     edm::Handle< L1GlobalTriggerReadoutRecord > L1GlobalTrigger;
787    
788     iEvent.getByLabel(L1gtReadout_, L1GlobalTrigger);
789     const TechnicalTriggerWord& technicalTriggerWordBeforeMask = L1GlobalTrigger->technicalTriggerWord();
790    
791     for (int i=0; i<64;i++)
792     {
793     jets_.l1TBit[i] = technicalTriggerWordBeforeMask.at(i);
794     }
795     jets_.nL1TBit = 64;
796    
797     int ntrigs = L1GlobalTrigger->decisionWord().size();
798     jets_.nL1ABit = ntrigs;
799    
800     for (int i=0; i != ntrigs; i++) {
801     bool accept = L1GlobalTrigger->decisionWord()[i];
802     //jets_.l1ABit[i] = (accept == true)? 1:0;
803     if(accept== true){
804     jets_.l1ABit[i] = 1;
805     }
806     else{
807     jets_.l1ABit[i] = 0;
808     }
809    
810     }
811     }
812    
813     //--------------------------------------------------------------------------------------------------
814     void HiInclusiveJetAnalyzer::fillHLTBits(const edm::Event &iEvent)
815     {
816     // Fill HLT trigger bits.
817     Handle<TriggerResults> triggerResultsHLT;
818     getProduct(hltResName_, triggerResultsHLT, iEvent);
819    
820     const TriggerResults *hltResults = triggerResultsHLT.product();
821     const TriggerNames & triggerNames = iEvent.triggerNames(*hltResults);
822    
823     jets_.nHLTBit = hltTrgNames_.size();
824    
825     for(size_t i=0;i<hltTrgNames_.size();i++){
826    
827     for(size_t j=0;j<triggerNames.size();++j) {
828    
829     if(triggerNames.triggerName(j) == hltTrgNames_[i]){
830    
831     //cout <<"hltTrgNames_(i) "<<hltTrgNames_[i]<<endl;
832     //cout <<"triggerName(j) "<<triggerNames.triggerName(j)<<endl;
833     //cout<<" result "<<triggerResultsHLT->accept(j)<<endl;
834     jets_.hltBit[i] = triggerResultsHLT->accept(j);
835     }
836    
837     }
838     }
839     }
840    
841     //--------------------------------------------------------------------------------------------------
842     template <typename TYPE>
843     inline void HiInclusiveJetAnalyzer::getProduct(const std::string name, edm::Handle<TYPE> &prod,
844     const edm::Event &event) const
845     {
846     // Try to access data collection from EDM file. We check if we really get just one
847     // product with the given name. If not we throw an exception.
848    
849     event.getByLabel(edm::InputTag(name),prod);
850     if (!prod.isValid())
851     throw edm::Exception(edm::errors::Configuration, "HiInclusiveJetAnalyzer::GetProduct()\n")
852     << "Collection with label '" << name << "' is not valid" << std::endl;
853     }
854    
855     //--------------------------------------------------------------------------------------------------
856     template <typename TYPE>
857     inline bool HiInclusiveJetAnalyzer::getProductSafe(const std::string name, edm::Handle<TYPE> &prod,
858     const edm::Event &event) const
859     {
860     // Try to safely access data collection from EDM file. We check if we really get just one
861     // product with the given name. If not, we return false.
862    
863     if (name.size()==0)
864     return false;
865    
866     try {
867     event.getByLabel(edm::InputTag(name),prod);
868     if (!prod.isValid())
869     return false;
870     } catch (...) {
871     return false;
872     }
873     return true;
874     }
875    
876    
877 mnguyen 1.8 int
878     HiInclusiveJetAnalyzer::getPFJetMuon(const pat::Jet& pfJet, const reco::PFCandidateCollection *pfCandidateColl)
879     {
880    
881     int pfMuonIndex = -1;
882     float ptMax = 0.;
883    
884    
885     for(unsigned icand=0;icand<pfCandidateColl->size(); icand++) {
886     const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
887    
888     int id = pfCandidate.particleId();
889     if(abs(id) != 3) continue;
890 yilmaz 1.1
891 mnguyen 1.8 if(reco::deltaR(pfJet,pfCandidate)>0.5) continue;
892    
893     double pt = pfCandidate.pt();
894     if(pt>ptMax){
895     ptMax = pt;
896     pfMuonIndex = (int) icand;
897     }
898     }
899    
900     return pfMuonIndex;
901    
902     }
903    
904    
905     double
906     HiInclusiveJetAnalyzer::getPtRel(const reco::PFCandidate lep, const pat::Jet& jet )
907     {
908    
909     float lj_x = jet.p4().px();
910     float lj_y = jet.p4().py();
911     float lj_z = jet.p4().pz();
912    
913     // absolute values squared
914     float lj2 = lj_x*lj_x+lj_y*lj_y+lj_z*lj_z;
915     float lep2 = lep.px()*lep.px()+lep.py()*lep.py()+lep.pz()*lep.pz();
916    
917     // projection vec(mu) to lepjet axis
918     float lepXlj = lep.px()*lj_x+lep.py()*lj_y+lep.pz()*lj_z;
919    
920     // absolute value squared and normalized
921     float pLrel2 = lepXlj*lepXlj/lj2;
922    
923     // lep2 = pTrel2 + pLrel2
924     float pTrel2 = lep2-pLrel2;
925    
926     return (pTrel2 > 0) ? std::sqrt(pTrel2) : 0.0;
927     }
928 mnguyen 1.11
929     // Recursive function, but this version gets called only the first time
930     void
931     HiInclusiveJetAnalyzer::saveDaughters(const reco::GenParticle &gen){
932    
933     for(unsigned i=0;i<gen.numberOfDaughters();i++){
934     const reco::Candidate & daughter = *gen.daughter(i);
935     double daughterPt = daughter.pt();
936     if(daughterPt<1.) continue;
937     double daughterEta = daughter.eta();
938     if(fabs(daughterEta)>3.) continue;
939     int daughterPdgId = daughter.pdgId();
940     int daughterStatus = daughter.status();
941     // Special case when b->b+string, both b and string contain all daughters, so only take the string
942     if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
943    
944     // cheesy way of finding strings which were already used
945     if(daughter.pdgId()==92){
946     for(unsigned ist=0;ist<usedStringPts.size();ist++){
947     if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
948     }
949     usedStringPts.push_back(daughter.pt());
950     }
951     jets_.bJetIndex[jets_.bMult] = jets_.nref;
952     jets_.bStatus[jets_.bMult] = daughterStatus;
953     jets_.bVx[jets_.bMult] = daughter.vx();
954     jets_.bVy[jets_.bMult] = daughter.vy();
955     jets_.bVz[jets_.bMult] = daughter.vz();
956     jets_.bPt[jets_.bMult] = daughterPt;
957     jets_.bEta[jets_.bMult] = daughterEta;
958     jets_.bPhi[jets_.bMult] = daughter.phi();
959     jets_.bPdg[jets_.bMult] = daughterPdgId;
960     jets_.bChg[jets_.bMult] = daughter.charge();
961     jets_.bMult++;
962     saveDaughters(daughter);
963     }
964     }
965    
966     // This version called for all subsequent calls
967     void
968     HiInclusiveJetAnalyzer::saveDaughters(const reco::Candidate &gen){
969    
970     for(unsigned i=0;i<gen.numberOfDaughters();i++){
971     const reco::Candidate & daughter = *gen.daughter(i);
972     double daughterPt = daughter.pt();
973     if(daughterPt<1.) continue;
974     double daughterEta = daughter.eta();
975     if(fabs(daughterEta)>3.) continue;
976     int daughterPdgId = daughter.pdgId();
977     int daughterStatus = daughter.status();
978     // Special case when b->b+string, both b and string contain all daughters, so only take the string
979     if(gen.pdgId()==daughterPdgId && gen.status()==3 && daughterStatus==2) continue;
980    
981     // cheesy way of finding strings which were already used
982     if(daughter.pdgId()==92){
983     for(unsigned ist=0;ist<usedStringPts.size();ist++){
984     if(fabs(daughter.pt() - usedStringPts[ist]) < 0.0001) return;
985     }
986     usedStringPts.push_back(daughter.pt());
987     }
988    
989     jets_.bJetIndex[jets_.bMult] = jets_.nref;
990     jets_.bStatus[jets_.bMult] = daughterStatus;
991     jets_.bVx[jets_.bMult] = daughter.vx();
992     jets_.bVy[jets_.bMult] = daughter.vy();
993     jets_.bVz[jets_.bMult] = daughter.vz();
994     jets_.bPt[jets_.bMult] = daughterPt;
995     jets_.bEta[jets_.bMult] = daughterEta;
996     jets_.bPhi[jets_.bMult] = daughter.phi();
997     jets_.bPdg[jets_.bMult] = daughterPdgId;
998     jets_.bChg[jets_.bMult] = daughter.charge();
999     jets_.bMult++;
1000     saveDaughters(daughter);
1001     }
1002     }
1003 yilmaz 1.1
1004     DEFINE_FWK_MODULE(HiInclusiveJetAnalyzer);