ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/JetAnalyzer.cc
Revision: 1.3
Committed: Mon Dec 15 19:08:52 2008 UTC (16 years, 4 months ago) by lethuill
Content type: text/plain
Branch: MAIN
CVS Tags: pat_2_1_12_01
Changes since 1.2: +217 -105 lines
Log Message:
MC infos from Silvano

File Contents

# User Rev Content
1 lethuill 1.2 #include "../interface/JetAnalyzer.h"
2 lethuill 1.1
3     using namespace std;
4     using namespace reco;
5     using namespace edm;
6    
7 lethuill 1.2 JetAnalyzer::JetAnalyzer(const edm::ParameterSet& producersNames):verbosity_(0),useMC_(false)
8 lethuill 1.1 {
9 lethuill 1.2 dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
10 lethuill 1.1 jetProducer_ = producersNames.getParameter<edm::InputTag>("jetProducer");
11     }
12    
13 lethuill 1.2 JetAnalyzer::JetAnalyzer(const edm::ParameterSet& producersNames, int verbosity):verbosity_(verbosity),useMC_(false)
14 lethuill 1.1 {
15 lethuill 1.2 dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
16 lethuill 1.1 jetProducer_ = producersNames.getParameter<edm::InputTag>("jetProducer");
17     }
18    
19 lethuill 1.2 JetAnalyzer::JetAnalyzer(const edm::ParameterSet& producersNames, const edm::ParameterSet& myConfig, int verbosity):verbosity_(verbosity)
20     {
21     dataType_ = producersNames.getUntrackedParameter<string>("dataType","unknown");
22     jetProducer_ = producersNames.getParameter<edm::InputTag>("jetProducer");
23 lethuill 1.3 useMC_ = myConfig.getUntrackedParameter<bool>("doJetMC");
24 lethuill 1.2 }
25    
26 lethuill 1.1 JetAnalyzer::~JetAnalyzer()
27     {
28     }
29    
30 lethuill 1.3 bool Rsortrule (std::pair <double,double> p1, std::pair <double,double> p2 )
31     {
32     return p1.second<p2.second;
33 lethuill 1.2 }
34    
35 lethuill 1.1 void JetAnalyzer::Process(const edm::Event& iEvent, TClonesArray* rootJets)
36     {
37    
38 lethuill 1.3 unsigned int nJets=0;
39    
40     // reco::CaloJet or reco::PFJet ?
41     std::string jetType = "BASIC";
42     if( jetProducer_.label()=="kt4CaloJets"
43     || jetProducer_.label()=="kt6CaloJets"
44     || jetProducer_.label()=="iterativeCone5CaloJets"
45     || jetProducer_.label()=="sisCone5CaloJets"
46     || jetProducer_.label()=="sisCone7CaloJets"
47     ) jetType="CALO";
48    
49     if( jetProducer_.label()=="kt4PFJets"
50     || jetProducer_.label()=="kt6PFJets"
51     || jetProducer_.label()=="iterativeCone5PFJets"
52     || jetProducer_.label()=="sisCone5PFJets"
53     || jetProducer_.label()=="sisCone7PFJets"
54     ) jetType="PF";
55    
56     edm::Handle < std::vector <reco::CaloJet> > recoCaloJets;
57     if( (dataType_=="RECO" || dataType_=="AOD") && jetType=="CALO" )
58     {
59     iEvent.getByLabel(jetProducer_, recoCaloJets);
60     nJets = recoCaloJets->size();
61     }
62    
63     edm::Handle < std::vector <reco::PFJet> > recoPFJets;
64     if( (dataType_=="RECO" || dataType_=="AOD") && jetType=="PF" )
65 lethuill 1.1 {
66 lethuill 1.3 iEvent.getByLabel(jetProducer_, recoPFJets);
67     nJets = recoPFJets->size();
68 lethuill 1.1 }
69 lethuill 1.3
70     edm::Handle < std::vector <pat::Jet> > patJets;
71     if( dataType_=="PAT" || dataType_=="PATAOD" )
72     {
73     iEvent.getByLabel(jetProducer_, patJets);
74     nJets = patJets->size();
75 lethuill 1.2 }
76    
77 lethuill 1.3 if(verbosity_>1) std::cout << " Number of jets = " << nJets << " Label: " << jetProducer_.label() << " Instance: " << jetProducer_.instance() << std::endl;
78    
79     for (unsigned int j=0; j<nJets; j++)
80     {
81     const reco::Jet* jet = 0;
82     if( (dataType_=="RECO" || dataType_=="AOD") && jetType=="CALO" ) jet = (const reco::Jet*) ( & ((*recoCaloJets)[j]) );
83     if( (dataType_=="RECO" || dataType_=="AOD") && jetType=="PF" ) jet = (const reco::Jet*) ( & ((*recoPFJets)[j]) );
84     if( dataType_=="PAT" || dataType_=="PATAOD" )
85     {
86     jet = (const reco::Jet*) ( & ((*patJets)[j]) );
87     if( (*patJets)[j].isCaloJet() ) jetType="CALO";
88     if( (*patJets)[j].isPFJet() ) jetType="PF";
89     }
90    
91     TRootJet localJet(
92     jet->px()
93     ,jet->py()
94     ,jet->pz()
95     ,jet->energy()
96     ,jet->vx()
97     ,jet->vy()
98     ,jet->vz()
99     ,abs(jet->pdgId())
100     ,jet->charge()
101     );
102    
103     localJet.setNConstituents(jet->nConstituents());
104     localJet.setN90(jet->nCarrying(0.9));
105     localJet.setN60(jet->nCarrying(0.6));
106     localJet.setJetArea(jet->jetArea());
107     localJet.setPileupEnergy(jet->pileup());
108     localJet.setMaxDistance(jet->maxDistance());
109     float totalEnergy = jet->etInAnnulus(0.,999.);
110     if(totalEnergy!=0)
111     {
112     localJet.setDR01EnergyFraction( jet->etInAnnulus(0.,0.1) / totalEnergy );
113     localJet.setDR02EnergyFraction( jet->etInAnnulus(0.,0.2) / totalEnergy );
114     localJet.setDR03EnergyFraction( jet->etInAnnulus(0.,0.3) / totalEnergy );
115     localJet.setDR04EnergyFraction( jet->etInAnnulus(0.,0.4) / totalEnergy );
116     localJet.setDR05EnergyFraction( jet->etInAnnulus(0.,0.5) / totalEnergy );
117     }
118    
119    
120     if( dataType_=="RECO" || dataType_=="AOD" )
121     {
122     // Some specific methods to reco::Jet
123     // Do association to genParticle ?
124    
125     if( jetType=="CALO" )
126     {
127     // Variables from reco::CaloJet
128     const reco::CaloJet *caloJet = dynamic_cast<const reco::CaloJet*>(&*jet);
129     localJet.setJetType(1);
130     localJet.setEcalEnergyFraction(caloJet->emEnergyFraction());
131     localJet.setHcalEnergyFraction(caloJet->energyFractionHadronic());
132     //std::vector<CaloTowerPtr> getCaloConstituents () const;
133     }
134    
135     if( jetType=="PF" )
136     {
137     // Variables from reco::PFJet
138     const reco::PFJet *pfJet = dynamic_cast<const reco::PFJet*>(&*jet);
139     localJet.setJetType(2);
140     localJet.setEcalEnergyFraction(pfJet->chargedEmEnergyFraction() + pfJet->neutralEmEnergyFraction());
141     localJet.setHcalEnergyFraction(pfJet->chargedHadronEnergyFraction() + pfJet->neutralHadronEnergyFraction());
142     if(pfJet->energy() != 0) localJet.setChargedEnergyFraction( (pfJet->chargedEmEnergy() + pfJet->chargedHadronEnergy() + pfJet->chargedMuEnergy() ) / pfJet->energy());
143     localJet.setChargedMultiplicity((int)pfJet->chargedMultiplicity()) ;
144     //std::vector <const reco::PFCandidate*> getPFConstituents () const;
145     }
146    
147     }
148    
149    
150     if( dataType_=="PATAOD" || dataType_=="PAT" )
151     {
152     // Some specific methods to pat::Jet
153     const pat::Jet *patJet = dynamic_cast<const pat::Jet*>(&*jet);
154    
155     // Variables from pat::Jet (Basic)
156     localJet.setBtag_trackCountingHighEff(patJet->bDiscriminator("trackCountingHighEffBJetTags"));
157     localJet.setBtag_trackCountingHighPur(patJet->bDiscriminator("trackCountingHighPurBJetTags"));
158     localJet.setBtag_jetProbability(patJet->bDiscriminator("jetProbabilityBJetTags"));
159     localJet.setBCorrection(patJet->correctionFactor(pat::Jet::bCorrection));
160     localJet.setCCorrection(patJet->correctionFactor(pat::Jet::cCorrection));
161     localJet.setUDSCorrection(patJet->correctionFactor(pat::Jet::udsCorrection));
162     localJet.setGCorrection(patJet->correctionFactor(pat::Jet::gCorrection));
163    
164     // Use associated tracks to calculate charged broadness of the jet
165     // FIXME - Check generalTracks collection is present
166     reco::TrackRefVector tracks = patJet->associatedTracks() ;
167     Float_t deltaR = 0.;
168     Float_t pTrackerTot = 0.;
169     Float_t pTrackerCone01 = 0.;
170     Float_t pTrackerCone02 = 0.;
171     Float_t pTrackerCone03 = 0.;
172     Float_t pTrackerCone04 = 0.;
173     Float_t pTrackerCone05 = 0.;
174     // TODO - Use std::map....
175     vector < pair < Float_t , Float_t > > tracksVPair ;
176     pair < Float_t , Float_t > tracksPair;
177    
178     for(unsigned int iTracks = 0; iTracks< patJet->associatedTracks().size() ; iTracks++)
179     {
180     deltaR = localJet.DeltaR(TLorentzVector(tracks[iTracks]->px(),tracks[iTracks]->py(),tracks[iTracks]->pz(),0));
181     pTrackerTot += tracks[iTracks]->p();
182     tracksPair.first = tracks[iTracks]->p();
183     tracksPair.second = deltaR;
184     tracksVPair.push_back(tracksPair);
185     if(deltaR < 0.5) pTrackerCone05 += tracks[iTracks]->p();
186     if(deltaR < 0.4) pTrackerCone04 += tracks[iTracks]->p();
187     if(deltaR < 0.3) pTrackerCone03 += tracks[iTracks]->p();
188     if(deltaR < 0.2) pTrackerCone02 += tracks[iTracks]->p();
189     if(deltaR < 0.1) pTrackerCone01 += tracks[iTracks]->p();
190     }
191     sort(tracksVPair.begin(), tracksVPair.end(), Rsortrule);
192     Float_t Rmin = 0;
193     Float_t pDummy = 0;
194     for(std::vector<std::pair< Float_t,Float_t > >::iterator i = tracksVPair.begin(); i!=tracksVPair.end(); i++)
195     {
196     pDummy+=(*i).first;
197     if (pDummy>0.75*(pTrackerTot))
198     {
199     Rmin = (*i).second;
200     break;
201     }
202     }
203     if (Rmin<1e-5) {Rmin=0.;}
204     localJet.setChargedBroadness(Rmin);
205     if(pTrackerTot !=0)
206     {
207     localJet.setChargedBroadnessDR01(pTrackerCone01/pTrackerTot);
208     localJet.setChargedBroadnessDR02(pTrackerCone02/pTrackerTot);
209     localJet.setChargedBroadnessDR03(pTrackerCone03/pTrackerTot);
210     localJet.setChargedBroadnessDR04(pTrackerCone04/pTrackerTot);
211     localJet.setChargedBroadnessDR05(pTrackerCone05/pTrackerTot);
212     }
213    
214    
215     if ( patJet->isCaloJet() ) {
216     // Variables from pat::Jet (CaloSpecific)
217     localJet.setJetType(1);
218     localJet.setEcalEnergyFraction(patJet->emEnergyFraction());
219     localJet.setHcalEnergyFraction(patJet->energyFractionHadronic());
220     localJet.setChargedMultiplicity(patJet->associatedTracks().size()) ;
221     //std::vector<CaloTowerPtr> getCaloConstituents () const;
222     }
223    
224     if(patJet->isPFJet()) {
225     // Variables from pat::Jet (PFSpecific)
226     localJet.setJetType(2);
227     localJet.setEcalEnergyFraction(patJet->chargedEmEnergyFraction() + patJet->neutralEmEnergyFraction());
228     localJet.setHcalEnergyFraction(patJet->chargedHadronEnergyFraction() + patJet->neutralHadronEnergyFraction());
229     if(patJet->energy() != 0) localJet.setChargedEnergyFraction( (patJet->chargedEmEnergy() + patJet->chargedHadronEnergy() + patJet->chargedMuEnergy() ) / patJet->energy());
230     localJet.setChargedMultiplicity((int)patJet->chargedMultiplicity()) ;
231     //std::vector <const reco::PFCandidate*> getPFConstituents () const;
232     }
233    
234     // Matched genParticle
235     if (useMC_)
236     {
237     if ( patJet->genParton() != NULL )
238     {
239     localJet.setMomentumMCParton(TLorentzVector(patJet->genParton()->px(),patJet->genParton()->py(),patJet->genParton()->pz(),patJet->genParton()->energy()));
240     localJet.setVertexMCParton(TVector3(patJet->genParton()->vx(),patJet->genParton()->vy(),patJet->genParton()->vz() ) );
241     localJet.setPdgIdMCParton(patJet->genParton()->pdgId());
242     }
243    
244     if ( patJet->genJet() != NULL )
245     {
246     localJet.setMomentumMCJet(TLorentzVector(patJet->genJet()->px(),patJet->genJet()->py(),patJet->genJet()->pz(),patJet->genJet()->energy()));
247     localJet.setVertexMCJet(TVector3(patJet->genJet()->vx(),patJet->genJet()->vy(),patJet->genJet()->vz() ));
248     localJet.setPdgIdMCJet(patJet->genJet()->pdgId());
249     }
250     }
251     }
252 lethuill 1.1
253 lethuill 1.3 new( (*rootJets)[j] ) TRootJet(localJet);
254     if(verbosity_>2) cout << " ["<< setw(3) << j << "] " << localJet << endl;
255     }
256 lethuill 1.1 }