ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Morgan/src/JetAnalyzer.cc
(Generate patch)

Comparing UserCode/Morgan/src/JetAnalyzer.cc (file contents):
Revision 1.9 by lethuill, Fri Apr 10 08:27:40 2009 UTC vs.
Revision 1.10 by lethuill, Wed Jun 10 11:17:06 2009 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines