ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerJPTJets.cc
Revision: 1.4
Committed: Wed Sep 22 08:39:34 2010 UTC (14 years, 7 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_032, Mit_031, Mit_025c_branch2, Mit_025c_branch1, Mit_030, Mit_029c, Mit_029b, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_025c_branch0, Mit_028, Mit_027a, Mit_027, Mit_026, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, HEAD
Branch point for: Mit_025c_branch
Changes since 1.3: +3 -3 lines
Log Message:
silence debug printout

File Contents

# User Rev Content
1 bendavid 1.4 // $Id: FillerJPTJets.cc,v 1.3 2010/09/19 23:47:33 bendavid Exp $
2 bendavid 1.1
3     #include "MitProd/TreeFiller/interface/FillerJPTJets.h"
4     #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
5     #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
6     #include "DataFormats/JetReco/interface/JPTJet.h"
7     #include "DataFormats/BTauReco/interface/JetTag.h"
8     #include "SimDataFormats/JetMatching/interface/JetFlavour.h"
9     #include "SimDataFormats/JetMatching/interface/JetFlavourMatching.h"
10     #include "SimDataFormats/JetMatching/interface/MatchedPartons.h"
11     #include "SimDataFormats/JetMatching/interface/JetMatchedPartons.h"
12     #include "JetMETCorrections/Objects/interface/JetCorrector.h"
13     #include "MitAna/DataTree/interface/Names.h"
14     #include "MitAna/DataTree/interface/CaloJetCol.h"
15     #include "MitAna/DataTree/interface/JPTJetCol.h"
16     #include "MitEdm/DataFormats/interface/RefToBaseToPtr.h"
17     #include "MitProd/ObjectService/interface/ObjectService.h"
18    
19     using namespace std;
20     using namespace edm;
21     using namespace mithep;
22    
23     //--------------------------------------------------------------------------------------------------
24     FillerJPTJets::FillerJPTJets(const ParameterSet &cfg, const char *name, bool active) :
25     BaseFiller(cfg,name,active),
26     flavorMatchingActive_(Conf().getUntrackedParameter<bool>("flavorMatchingActive",true)),
27     bTaggingActive_(Conf().getUntrackedParameter<bool>("bTaggingActive",true)),
28     jetToVertexActive_(Conf().getUntrackedParameter<bool>("jetToVertexActive",true)),
29     jetCorrectionsActive_(Conf().getUntrackedParameter<bool>("jetCorrectionsActive",true)),
30     edmName_(Conf().getUntrackedParameter<string>("edmName","recoJPTJets:iterativeCone5JPTJets")),
31     mitName_(Conf().getUntrackedParameter<string>("mitName","ItrCone5JPTJets")),
32     jetToVertexAlphaName_(Conf().getUntrackedParameter<string>
33     ("jetToVertexAlphaName","jetToVertexAlpha")),
34     jetToVertexBetaName_(Conf().getUntrackedParameter<string>
35     ("jetToVertexBetaName","jetToVertexBetaName")),
36     L2JetCorrectorName_(Conf().getUntrackedParameter<string>
37     ("L2JetCorrectorName","L2JetCorrectorName")),
38     L3JetCorrectorName_(Conf().getUntrackedParameter<string>
39     ("L3JetCorrectorName","L3JetCorrectorName")),
40     flavorMatchingByReferenceName_(Conf().getUntrackedParameter<string>
41     ("flavorMatchingByReferenceName","srcByReference")),
42     flavorMatchingDefinition_(Conf().getUntrackedParameter<string>
43     ("flavorMatchingDefinition","Algorithmic")),
44     jetProbabilityBJetTagsName_(Conf().getUntrackedParameter<string>
45     ("JetProbabilityBJetTagsName","jetProbabilityBJetTags")),
46     jetBProbabilityBJetTagsName_(Conf().getUntrackedParameter<string>
47     ("JetBProbabilityBJetTagsName","jetBProbabilityBJetTags")),
48     simpleSecondaryVertexBJetTagsName_(Conf().getUntrackedParameter<string>
49     ("SimpleSecondaryVertexBJetTagsName","simpleSecondaryVertexBJetTags")),
50     combinedSecondaryVertexBJetTagsName_(Conf().getUntrackedParameter<string>
51     ("CombinedSecondaryVertexBJetTagsName","combinedSecondaryVertexBJetTags")),
52     combinedSecondaryVertexMVABJetTagsName_(Conf().getUntrackedParameter<string>
53     ("CombinedSecondaryVertexMVABJetTagsName","combinedSecondaryVertexMVABJetTags")),
54     trackCountingHighEffBJetTagsName_(Conf().getUntrackedParameter<string>
55     ("TrackCountingHighEffBJetTagsName","trackCountingHighEffBJetTags")),
56     trackCountingHighPurBJetTagsName_(Conf().getUntrackedParameter<string>
57     ("TrackCountingHighPurBJetTagsName","trackCountingHighPurBJetTags")),
58     softMuonBJetTagsName_(Conf().getUntrackedParameter<string>
59     ("SoftMuonBJetTagsName","softMuonBJetTags")),
60     softMuonByIP3dBJetTagsName_(Conf().getUntrackedParameter<string>
61     ("SoftMuonByIP3dBJetTagsName","softMuonByIP3dBJetTags")),
62     softMuonByPtBJetTagsName_(Conf().getUntrackedParameter<string>
63     ("SoftMuonByPtBJetTagsName","softMuonByPtBJetTags")),
64     softElectronByIP3dBJetTagsName_(Conf().getUntrackedParameter<string>
65     ("SoftElectronByIP3dBJetTagsName","softElectronByIP3dBJetTags")),
66     softElectronByPtBJetTagsName_(Conf().getUntrackedParameter<string>
67     ("SoftElectronByPtBJetTagsName","softElectronByPtBJetTags")),
68     caloJetMapName_(Conf().getUntrackedParameter<string>("caloJetCandMapName","caloJetMapName")),
69     jetMapName_(Conf().getUntrackedParameter<string>("jetMapName","JPTJetMap")),
70     caloJetMap_(0),
71     jetMap_(new mithep::JPTJetMap),
72     jets_(new mithep::JPTJetArr(16))
73     {
74     // Constructor.
75     }
76    
77     //--------------------------------------------------------------------------------------------------
78     FillerJPTJets::~FillerJPTJets()
79     {
80     // Destructor.
81    
82     delete jets_;
83     delete jetMap_;
84     }
85    
86     //--------------------------------------------------------------------------------------------------
87     void FillerJPTJets::BookDataBlock(TreeWriter &tws)
88     {
89     // Add jets branch to tree.
90    
91     tws.AddBranch(mitName_,&jets_);
92     OS()->add<mithep::JPTJetArr>(jets_,mitName_);
93    
94     if (!caloJetMapName_.empty()) {
95     caloJetMap_ = OS()->get<CaloJetMap>(caloJetMapName_);
96     if (caloJetMap_)
97     AddBranchDep(mitName_,caloJetMap_->GetBrName());
98     }
99     if (!jetMapName_.empty()) {
100     jetMap_->SetBrName(mitName_);
101     OS()->add<JPTJetMap>(jetMap_,jetMapName_);
102     }
103     }
104    
105     //--------------------------------------------------------------------------------------------------
106     void FillerJPTJets::FillDataBlock(const edm::Event &event,
107     const edm::EventSetup &setup)
108     {
109     // Fill jets from edm collection into our collection.
110    
111     jets_->Delete();
112     jetMap_->Reset();
113 bendavid 1.3
114     //ugly hack to get collection from original reco if present (needed to handle things
115     //automatically between 35x and post 35x samples in the correct way)
116    
117     Handle<reco::PFJetCollection> hPFJetProduct;
118     event.getByLabel("ak5PFJets",hPFJetProduct);
119 bendavid 1.1
120     // handle for the Jet Collection
121     Handle<reco::JPTJetCollection> hJetProduct;
122 bendavid 1.3
123     if (hPFJetProduct.isValid()) {
124 bendavid 1.4 //printf("Getting jpt from original reco\n");
125 bendavid 1.3 InputTag jpttag(edmName_,"",hPFJetProduct.provenance()->processName());
126     event.getByLabel(jpttag,hJetProduct);
127     }
128    
129     if (!hPFJetProduct.isValid() || !hJetProduct.isValid()) {
130 bendavid 1.4 //printf("falling back to newly produced jpt collection\n");
131 bendavid 1.3 event.getByLabel(edmName_,hJetProduct);
132     }
133    
134    
135 bendavid 1.1 // handles for jet flavour matching
136     Handle<reco::JetMatchedPartonsCollection> hPartonMatchingProduct;
137     if (flavorMatchingActive_)
138     GetProduct(flavorMatchingByReferenceName_, hPartonMatchingProduct, event);
139    
140     Handle<reco::JetTagCollection> hJetProbabilityBJetTags;
141     Handle<reco::JetTagCollection> hJetBProbabilityBJetTags;
142     Handle<reco::JetTagCollection> hSimpleSecondaryVertexBJetTags;
143     Handle<reco::JetTagCollection> hCombinedSecondaryVertexBJetTags;
144     Handle<reco::JetTagCollection> hCombinedSecondaryVertexMVABJetTags;
145     Handle<reco::JetTagCollection> hTrackCountingHighEffBJetTags;
146     Handle<reco::JetTagCollection> hTrackCountingHighPurBJetTags;
147     Handle<reco::JetTagCollection> hSoftMuonBJetTags;
148     Handle<reco::JetTagCollection> hSoftMuonByIP3dBJetTags;
149     Handle<reco::JetTagCollection> hSoftMuonByPtBJetTags;
150     Handle<reco::JetTagCollection> hSoftElectronByIP3dBJetTags;
151     Handle<reco::JetTagCollection> hSoftElectronByPtBJetTags;
152    
153     if (bTaggingActive_) {
154     GetProduct(jetProbabilityBJetTagsName_, hJetProbabilityBJetTags, event);
155     GetProduct(jetBProbabilityBJetTagsName_, hJetBProbabilityBJetTags, event);
156     GetProduct(simpleSecondaryVertexBJetTagsName_, hSimpleSecondaryVertexBJetTags, event);
157     GetProduct(combinedSecondaryVertexBJetTagsName_, hCombinedSecondaryVertexBJetTags, event);
158     GetProduct(combinedSecondaryVertexMVABJetTagsName_, hCombinedSecondaryVertexMVABJetTags, event);
159     GetProduct(trackCountingHighEffBJetTagsName_, hTrackCountingHighEffBJetTags, event);
160     GetProduct(trackCountingHighPurBJetTagsName_, hTrackCountingHighPurBJetTags, event);
161     GetProduct(softMuonBJetTagsName_, hSoftMuonBJetTags, event);
162     GetProduct(softMuonByIP3dBJetTagsName_, hSoftMuonByIP3dBJetTags, event);
163     GetProduct(softMuonByPtBJetTagsName_, hSoftMuonByPtBJetTags, event);
164     GetProduct(softElectronByIP3dBJetTagsName_, hSoftElectronByIP3dBJetTags, event);
165     GetProduct(softElectronByPtBJetTagsName_, hSoftElectronByPtBJetTags, event);
166     }
167    
168     const reco::JPTJetCollection inJets = *(hJetProduct.product());
169    
170     //Handles to Jet to Vertex Association
171     Handle<std::vector<double> > JV_alpha;
172     Handle<std::vector<double> > JV_beta;
173     std::vector<double>::const_iterator it_jv_alpha;
174     std::vector<double>::const_iterator it_jv_beta;
175    
176     if (jetToVertexActive_) {
177     GetProduct(jetToVertexAlphaName_, JV_alpha, event);
178     GetProduct(jetToVertexBetaName_, JV_beta, event);
179     it_jv_alpha = JV_alpha->begin();
180     it_jv_beta = JV_beta->begin();
181     }
182    
183     //Define Jet Correction Services
184     const JetCorrector* correctorL2 = 0;
185     const JetCorrector* correctorL3 = 0;
186     if (jetCorrectionsActive_) {
187     correctorL2 = JetCorrector::getJetCorrector (L2JetCorrectorName_,setup);
188     correctorL3 = JetCorrector::getJetCorrector (L3JetCorrectorName_,setup);
189     }
190    
191     // loop through all jets
192     for (reco::JPTJetCollection::const_iterator inJet = inJets.begin();
193     inJet != inJets.end(); ++inJet) {
194    
195     reco::JPTJetRef jetRef(hJetProduct, inJet - inJets.begin());
196     reco::JetBaseRef jetBaseRef(jetRef);
197    
198     mithep::JPTJet *jet = jets_->Allocate();
199     new (jet) mithep::JPTJet(inJet->p4().x(),
200     inJet->p4().y(),
201     inJet->p4().z(),
202     inJet->p4().e());
203    
204     //add to map
205     edm::Ptr<reco::JPTJet> thePtr(hJetProduct, inJet - inJets.begin());
206     jetMap_->Add(thePtr, jet);
207    
208     //fill jet moments
209     jet->SetSigmaEta(TMath::Sqrt(inJet->etaetaMoment()));
210     jet->SetSigmaPhi(TMath::Sqrt(inJet->phiphiMoment()));
211    
212    
213     //fill jptjet-specific quantities
214     jet->SetZSPCor(inJet->getZSPCor());
215     jet->SetChargedHadronEnergy(inJet->chargedHadronEnergy());
216     jet->SetNeutralHadronEnergy(inJet->neutralHadronEnergy());
217     jet->SetChargedEmEnergy(inJet->chargedHadronEnergy());
218     jet->SetNeutralEmEnergy(inJet->neutralEmEnergy());
219     jet->SetResponseOfChargedWithEff(inJet->getSpecific().mResponseOfChargedWithEff);
220     jet->SetResponseOfChargedWithoutEff(inJet->getSpecific().mResponseOfChargedWithoutEff);
221     jet->SetSumPtOfChargedWithEff(inJet->getSpecific().mSumPtOfChargedWithEff);
222     jet->SetSumPtOfChargedWithoutEff(inJet->getSpecific().mSumPtOfChargedWithoutEff);
223     jet->SetSumEnergyOfChargedWithEff(inJet->getSpecific().mSumEnergyOfChargedWithEff);
224     jet->SetSumEnergyOfChargedWithoutEff(inJet->getSpecific().mSumEnergyOfChargedWithoutEff);
225     jet->SetR2momtr(inJet->getSpecific().R2momtr);
226     jet->SetEta2momtr(inJet->getSpecific().Eta2momtr);
227     jet->SetPhi2momtr(inJet->getSpecific().Phi2momtr);
228     jet->SetPout(inJet->getSpecific().Pout);
229     jet->SetZch(inJet->getSpecific().Zch);
230     jet->SetChargedMultiplicity(inJet->chargedMultiplicity());
231     jet->SetMuonMultiplicity(inJet->muonMultiplicity());
232     jet->SetElectronMultiplicity(inJet->elecMultiplicity());
233    
234     if (jetToVertexActive_) {
235     //compute alpha and beta parameter for jets
236     jet->SetAlpha((*it_jv_alpha));
237     jet->SetBeta((*it_jv_beta));
238     }
239    
240     //Jet Corrections
241     if (jetCorrectionsActive_) {
242     double L2Scale = correctorL2->correction(inJet->p4());
243     double L3Scale = correctorL3->correction(inJet->p4()*L2Scale);
244     jet->SetL2RelativeCorrectionScale(L2Scale);
245     jet->SetL3AbsoluteCorrectionScale(L3Scale);
246     jet->EnableCorrection(mithep::JPTJet::L2);
247     jet->EnableCorrection(mithep::JPTJet::L3);
248     }
249    
250     if (bTaggingActive_) {
251     jet->SetJetProbabilityBJetTagsDisc((*(hJetProbabilityBJetTags.product()))[jetBaseRef]);
252     jet->SetJetBProbabilityBJetTagsDisc((*(hJetBProbabilityBJetTags.product()))[jetBaseRef]);
253     jet->SetSimpleSecondaryVertexBJetTagsDisc(
254     (*(hSimpleSecondaryVertexBJetTags.product()))[jetBaseRef]);
255     jet->SetCombinedSecondaryVertexBJetTagsDisc(
256     (*(hCombinedSecondaryVertexBJetTags.product()))[jetBaseRef]);
257     jet->SetCombinedSecondaryVertexMVABJetTagsDisc(
258     (*(hCombinedSecondaryVertexMVABJetTags.product()))[jetBaseRef]);
259     jet->SetTrackCountingHighEffBJetTagsDisc(
260     (*(hTrackCountingHighEffBJetTags.product()))[jetBaseRef]);
261     jet->SetTrackCountingHighPurBJetTagsDisc(
262     (*(hTrackCountingHighPurBJetTags.product()))[jetBaseRef]);
263     jet->SetSoftMuonBJetTagsDisc((*(hSoftMuonBJetTags.product()))[jetBaseRef]);
264     jet->SetSoftMuonByIP3dBJetTagsDisc((*(hSoftMuonByIP3dBJetTags.product()))[jetBaseRef]);
265     jet->SetSoftMuonByPtBJetTagsDisc((*(hSoftMuonByPtBJetTags.product()))[jetBaseRef]);
266     jet->SetSoftElectronByIP3dBJetTagsDisc((*(hSoftElectronByIP3dBJetTags.product()))[jetBaseRef]);
267     jet->SetSoftElectronByPtBJetTagsDisc((*(hSoftElectronByPtBJetTags.product()))[jetBaseRef]);
268     }
269    
270     // get the Monte Carlo flavour matching information
271     if (flavorMatchingActive_) {
272     unsigned int iJet = inJet - inJets.begin();
273     const reco::JetMatchedPartonsCollection *matchedPartons = hPartonMatchingProduct.product();
274     reco::MatchedPartons matchedParton = (*matchedPartons)[matchedPartons->key(iJet)];
275     int flavorPhysDef = (matchedParton.physicsDefinitionParton().isNonnull())?
276     matchedParton.physicsDefinitionParton()->pdgId():0;
277     int flavorAlgDef = (matchedParton.algoDefinitionParton().isNonnull())?
278     matchedParton.algoDefinitionParton()->pdgId():0;
279    
280     if (flavorMatchingDefinition_ == "Algorithmic") {
281     jet->SetMatchedMCFlavor(flavorAlgDef);
282     } else if(flavorMatchingDefinition_ == "Physics") {
283     jet->SetMatchedMCFlavor(flavorPhysDef);
284     } else {
285     jet->SetMatchedMCFlavor(0);
286     }
287     }
288    
289     // fill ref to original jet
290 bendavid 1.2 if (caloJetMap_ && inJet->getCaloJetRef().isNonnull()) {
291 bendavid 1.1 const edm::Ptr<reco::CaloJet> caloJetPtr(mitedm::refToBaseToPtr(inJet->getCaloJetRef()));
292     const mithep::Jet *originalJet = caloJetMap_->GetMit(caloJetPtr);
293     jet->SetOriginalJet(originalJet);
294     }
295     if (jetToVertexActive_) {
296     it_jv_alpha++;
297     it_jv_beta++;
298     }
299     }
300     jets_->Trim();
301     }