ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/MyVirtualJetProducer.cc
Revision: 1.7
Committed: Wed Sep 7 14:07:14 2011 UTC (13 years, 7 months ago) by frankma
Content type: text/plain
Branch: MAIN
CVS Tags: hi44X_02, hi413_03, hi441_1, hi441_0, hi413_11, hi413_10, hi413_09, hi413_08, hi413_07, hi413_06, hi413_05, hi413_04, hi413_02, hi39X_01, tag_d20110915
Changes since 1.6: +1 -1 lines
Log Message:
update to 44X

File Contents

# User Rev Content
1 yilmaz 1.1 ////////////////////////////////////////////////////////////////////////////////
2     //
3     // MyVirtualJetProducer
4     // ------------------
5     //
6     // 04/21/2009 Philipp Schieferdecker <philipp.schieferdecker@cern.ch>
7     ////////////////////////////////////////////////////////////////////////////////
8    
9     #include "MyVirtualJetProducer.h"
10     #include "RecoJets/JetProducers/interface/JetSpecific.h"
11    
12     #include "FWCore/Framework/interface/Event.h"
13     #include "FWCore/Framework/interface/EventSetup.h"
14     #include "FWCore/Framework/interface/ESHandle.h"
15 frankma 1.7 #include "FWCore/Utilities/interface/Exception.h"
16 yilmaz 1.1 #include "FWCore/MessageLogger/interface/MessageLogger.h"
17     #include "FWCore/Framework/interface/MakerMacros.h"
18    
19     #include "DataFormats/Common/interface/View.h"
20     #include "DataFormats/Common/interface/Handle.h"
21     #include "DataFormats/VertexReco/interface/Vertex.h"
22     #include "DataFormats/VertexReco/interface/VertexFwd.h"
23     #include "DataFormats/JetReco/interface/CaloJetCollection.h"
24     #include "DataFormats/JetReco/interface/GenJetCollection.h"
25     #include "DataFormats/JetReco/interface/PFJetCollection.h"
26     #include "DataFormats/JetReco/interface/BasicJetCollection.h"
27     #include "DataFormats/JetReco/interface/TrackJetCollection.h"
28     #include "DataFormats/Candidate/interface/CandidateFwd.h"
29     #include "DataFormats/Candidate/interface/LeafCandidate.h"
30 yilmaz 1.2 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
31     #include "Geometry/Records/interface/CaloGeometryRecord.h"
32    
33     #include "DataFormats/Math/interface/Vector3D.h"
34 yilmaz 1.1
35     #include "fastjet/SISConePlugin.hh"
36     #include "fastjet/CMSIterativeConePlugin.hh"
37     #include "fastjet/ATLASConePlugin.hh"
38     #include "fastjet/CDFMidPointPlugin.hh"
39    
40     #include <iostream>
41     #include <memory>
42     #include <algorithm>
43     #include <limits>
44     #include <cmath>
45    
46     using namespace std;
47    
48    
49     namespace reco {
50     namespace helper {
51     struct GreaterByPtPseudoJet {
52     bool operator()( const fastjet::PseudoJet & t1, const fastjet::PseudoJet & t2 ) const {
53     return t1.perp() > t2.perp();
54     }
55     };
56    
57     }
58     }
59    
60     //______________________________________________________________________________
61     const char *MyVirtualJetProducer::JetType::names[] = {
62     "BasicJet","GenJet","CaloJet","PFJet","TrackJet"
63     };
64    
65    
66     //______________________________________________________________________________
67     MyVirtualJetProducer::JetType::Type
68     MyVirtualJetProducer::JetType::byName(const string &name)
69     {
70     const char **pos = std::find(names, names + LastJetType, name);
71     if (pos == names + LastJetType) {
72     std::string errorMessage="Requested jetType not supported: "+name+"\n";
73     throw cms::Exception("Configuration",errorMessage);
74     }
75     return (Type)(pos-names);
76     }
77    
78    
79     void MyVirtualJetProducer::makeProduces( std::string alias, std::string tag )
80     {
81     if (makeCaloJet(jetTypeE)) {
82     produces<reco::CaloJetCollection>(tag).setBranchAlias(alias);
83     }
84     else if (makePFJet(jetTypeE)) {
85     produces<reco::PFJetCollection>(tag).setBranchAlias(alias);
86     }
87     else if (makeGenJet(jetTypeE)) {
88     produces<reco::GenJetCollection>(tag).setBranchAlias(alias);
89     }
90     else if (makeTrackJet(jetTypeE)) {
91     produces<reco::TrackJetCollection>(tag).setBranchAlias(alias);
92     }
93     else if (makeBasicJet(jetTypeE)) {
94     produces<reco::BasicJetCollection>(tag).setBranchAlias(alias);
95     }
96     }
97    
98     ////////////////////////////////////////////////////////////////////////////////
99     // construction / destruction
100     ////////////////////////////////////////////////////////////////////////////////
101    
102     //______________________________________________________________________________
103     MyVirtualJetProducer::MyVirtualJetProducer(const edm::ParameterSet& iConfig)
104     : moduleLabel_ (iConfig.getParameter<string> ("@module_label"))
105     , src_ (iConfig.getParameter<edm::InputTag>("src"))
106     , srcPVs_ (iConfig.getParameter<edm::InputTag>("srcPVs"))
107     , jetType_ (iConfig.getParameter<string> ("jetType"))
108     , jetAlgorithm_ (iConfig.getParameter<string> ("jetAlgorithm"))
109     , rParam_ (iConfig.getParameter<double> ("rParam"))
110     , inputEtMin_ (iConfig.getParameter<double> ("inputEtMin"))
111     , inputEMin_ (iConfig.getParameter<double> ("inputEMin"))
112     , jetPtMin_ (iConfig.getParameter<double> ("jetPtMin"))
113     , doPVCorrection_(iConfig.getParameter<bool> ("doPVCorrection"))
114     , restrictInputs_(false)
115     , maxInputs_(99999999)
116     , doAreaFastjet_ (iConfig.getParameter<bool> ("doAreaFastjet"))
117     , doRhoFastjet_ (iConfig.getParameter<bool> ("doRhoFastjet"))
118     , doPUOffsetCorr_(iConfig.getParameter<bool> ("doPUOffsetCorr"))
119     , maxBadEcalCells_ (iConfig.getParameter<unsigned>("maxBadEcalCells"))
120     , maxRecoveredEcalCells_ (iConfig.getParameter<unsigned>("maxRecoveredEcalCells"))
121     , maxProblematicEcalCells_(iConfig.getParameter<unsigned>("maxProblematicEcalCells"))
122     , maxBadHcalCells_ (iConfig.getParameter<unsigned>("maxBadHcalCells"))
123     , maxRecoveredHcalCells_ (iConfig.getParameter<unsigned>("maxRecoveredHcalCells"))
124     , maxProblematicHcalCells_(iConfig.getParameter<unsigned>("maxProblematicHcalCells"))
125     , jetCollInstanceName_ ("")
126     {
127 yilmaz 1.2 // ntuple = new TNtuple("nt","debug","ieta:eta:iphi:phi");
128    
129 yilmaz 1.1 //
130     // additional parameters to think about:
131     // - overlap threshold (set to 0.75 for the time being)
132     // - p parameter for generalized kT (set to -2 for the time being)
133     // - fastjet PU subtraction parameters (not yet considered)
134     //
135     if (jetAlgorithm_=="SISCone") {
136     fjPlugin_ = PluginPtr( new fastjet::SISConePlugin(rParam_,0.75,0,0.0,false,
137     fastjet::SISConePlugin::SM_pttilde) );
138     fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(&*fjPlugin_) );
139     }
140     else if (jetAlgorithm_=="IterativeCone") {
141     fjPlugin_ = PluginPtr(new fastjet::CMSIterativeConePlugin(rParam_,1.0));
142     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(&*fjPlugin_));
143     }
144     else if (jetAlgorithm_=="CDFMidPoint") {
145     fjPlugin_ = PluginPtr(new fastjet::CDFMidPointPlugin(rParam_,0.75));
146     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(&*fjPlugin_));
147     }
148     else if (jetAlgorithm_=="ATLASCone") {
149     fjPlugin_ = PluginPtr(new fastjet::ATLASConePlugin(rParam_));
150     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(&*fjPlugin_));
151     }
152     else if (jetAlgorithm_=="Kt")
153     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(fastjet::kt_algorithm,rParam_));
154     else if (jetAlgorithm_=="CambridgeAachen")
155     fjJetDefinition_= JetDefPtr(new fastjet::JetDefinition(fastjet::cambridge_algorithm,
156     rParam_) );
157     else if (jetAlgorithm_=="AntiKt")
158     fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::antikt_algorithm,rParam_) );
159     else if (jetAlgorithm_=="GeneralizedKt")
160     fjJetDefinition_= JetDefPtr( new fastjet::JetDefinition(fastjet::genkt_algorithm,
161     rParam_,-2) );
162     else
163     throw cms::Exception("Invalid jetAlgorithm")
164     <<"Jet algorithm for MyVirtualJetProducer is invalid, Abort!\n";
165    
166     jetTypeE=JetType::byName(jetType_);
167    
168     if ( iConfig.exists("jetCollInstanceName") ) {
169     jetCollInstanceName_ = iConfig.getParameter<string>("jetCollInstanceName");
170     }
171    
172     if ( doPUOffsetCorr_ ) {
173 yilmaz 1.3 if ( jetTypeE != JetType::CaloJet ) {
174     throw cms::Exception("InvalidInput") << "Can only offset correct jets of type CaloJet";
175     }
176 yilmaz 1.4
177 yilmaz 1.3 puSubtractorName_ = iConfig.getParameter<string> ("subtractorName");
178 yilmaz 1.5
179     cout<<"puSubtractorName_ is : "<<puSubtractorName_.data()<<endl;
180 yilmaz 1.3 if(puSubtractorName_.empty()){
181     edm::LogWarning("VirtualJetProducer") << "Pile Up correction on; however, pile up type is not specified. Using default... \n";
182 yilmaz 1.5 cout<<"Pile Up correction on; however, pile up type is not specified. Using default A.K... \n";
183 yilmaz 1.3 subtractor_ = boost::shared_ptr<PileUpSubtractor>(new PileUpSubtractor(iConfig));
184 yilmaz 1.5
185 yilmaz 1.3 }else{
186 yilmaz 1.5 cout<<"getting subtractor "<<endl;
187 yilmaz 1.3 subtractor_ = boost::shared_ptr<PileUpSubtractor>(PileUpSubtractorFactory::get()->create( puSubtractorName_, iConfig));
188     }
189 yilmaz 1.1 }
190    
191     // do fasjet area / rho calcluation? => accept corresponding parameters
192     if ( doAreaFastjet_ || doRhoFastjet_ ) {
193     // default Ghost_EtaMax should be 5
194     double ghostEtaMax = iConfig.getParameter<double>("Ghost_EtaMax");
195     // default Active_Area_Repeats 1
196     int activeAreaRepeats = iConfig.getParameter<int> ("Active_Area_Repeats");
197     // default GhostArea 0.01
198     double ghostArea = iConfig.getParameter<double> ("GhostArea");
199     fjActiveArea_ = ActiveAreaSpecPtr(new fastjet::ActiveAreaSpec(ghostEtaMax,
200     activeAreaRepeats,
201     ghostArea));
202     fjRangeDef_ = RangeDefPtr( new fastjet::RangeDefinition(ghostEtaMax) );
203 yilmaz 1.4
204     }
205    
206     if ( doRhoFastjet_ ) {
207     doFastJetNonUniform_ = iConfig.getParameter<bool> ("doFastJetNonUniform");
208     if(doFastJetNonUniform_){
209     puCenters_ = iConfig.getParameter<std::vector<double> >("puCenters");
210     puWidth_ = iConfig.getParameter<double>("puWidth");
211     produces<std::vector<double> >("rhos");
212     produces<std::vector<double> >("sigmas");
213     }else{
214     produces<double>("rho");
215     produces<double>("sigma");
216     }
217    
218     }
219 yilmaz 1.1
220     // restrict inputs to first "maxInputs" towers?
221     if ( iConfig.exists("restrictInputs") ) {
222     restrictInputs_ = iConfig.getParameter<bool>("restrictInputs");
223     maxInputs_ = iConfig.getParameter<unsigned int>("maxInputs");
224     }
225    
226    
227     string alias=iConfig.getUntrackedParameter<string>("alias",moduleLabel_);
228    
229     // make the "produces" statements
230     makeProduces( alias, jetCollInstanceName_ );
231    
232     if(doRhoFastjet_){
233     produces<double>("rho");
234     produces<double>("sigma");
235     }
236     }
237    
238    
239     //______________________________________________________________________________
240     MyVirtualJetProducer::~MyVirtualJetProducer()
241     {
242     }
243    
244    
245     ////////////////////////////////////////////////////////////////////////////////
246     // implementation of member functions
247     ////////////////////////////////////////////////////////////////////////////////
248    
249     //______________________________________________________________________________
250     void MyVirtualJetProducer::produce(edm::Event& iEvent,const edm::EventSetup& iSetup)
251     {
252 yilmaz 1.2
253     if(!geo){
254     edm::ESHandle<CaloGeometry> pGeo;
255     iSetup.get<CaloGeometryRecord>().get(pGeo);
256     geo = pGeo.product();
257     }
258    
259 yilmaz 1.1 LogDebug("MyVirtualJetProducer") << "Entered produce\n";
260     //determine signal vertex
261     vertex_=reco::Jet::Point(0,0,0);
262     if (makeCaloJet(jetTypeE)&&doPVCorrection_) {
263     LogDebug("MyVirtualJetProducer") << "Adding PV info\n";
264     edm::Handle<reco::VertexCollection> pvCollection;
265     iEvent.getByLabel(srcPVs_,pvCollection);
266     if (pvCollection->size()>0) vertex_=pvCollection->begin()->position();
267     }
268    
269     // For Pileup subtraction using offset correction:
270     // set up geometry map
271     if ( doPUOffsetCorr_ ) {
272     subtractor_->setupGeometryMap(iEvent, iSetup);
273     }
274    
275     // clear data
276     LogDebug("MyVirtualJetProducer") << "Clear data\n";
277     fjInputs_.clear();
278     fjJets_.clear();
279     inputs_.clear();
280    
281     // get inputs and convert them to the fastjet format (fastjet::PeudoJet)
282     edm::Handle<reco::CandidateView> inputsHandle;
283     iEvent.getByLabel(src_,inputsHandle);
284     for (size_t i = 0; i < inputsHandle->size(); ++i) {
285     inputs_.push_back(inputsHandle->ptrAt(i));
286     }
287     LogDebug("MyVirtualJetProducer") << "Got inputs\n";
288    
289     // Convert candidates to fastjet::PseudoJets.
290     // Also correct to Primary Vertex. Will modify fjInputs_
291     // and use inputs_
292     fjInputs_.reserve(inputs_.size());
293     inputTowers();
294     LogDebug("MyVirtualJetProducer") << "Inputted towers\n";
295    
296     // For Pileup subtraction using offset correction:
297     // Subtract pedestal.
298     if ( doPUOffsetCorr_ ) {
299 yilmaz 1.4 subtractor_->reset(inputs_,fjInputs_,fjJets_);
300 yilmaz 1.1 subtractor_->calculatePedestal(fjInputs_);
301     subtractor_->subtractPedestal(fjInputs_);
302     LogDebug("MyVirtualJetProducer") << "Subtracted pedestal\n";
303     }
304    
305     // Run algorithm. Will modify fjJets_ and allocate fjClusterSeq_.
306     // This will use fjInputs_
307     runAlgorithm( iEvent, iSetup );
308     if ( doPUOffsetCorr_ ) {
309     subtractor_->setAlgorithm(fjClusterSeq_);
310     }
311    
312     LogDebug("MyVirtualJetProducer") << "Ran algorithm\n";
313    
314     // For Pileup subtraction using offset correction:
315     // Now we find jets and need to recalculate their energy,
316     // mark towers participated in jet,
317     // remove occupied towers from the list and recalculate mean and sigma
318     // put the initial towers collection to the jet,
319     // and subtract from initial towers in jet recalculated mean and sigma of towers
320     if ( doPUOffsetCorr_ ) {
321     vector<fastjet::PseudoJet> orphanInput;
322     subtractor_->calculateOrphanInput(orphanInput);
323     subtractor_->calculatePedestal(orphanInput);
324     subtractor_->offsetCorrectJets();
325     }
326    
327     // Write the output jets.
328     // This will (by default) call the member function template
329     // "writeJets", but can be overridden.
330     // this will use inputs_
331     output( iEvent, iSetup );
332     LogDebug("MyVirtualJetProducer") << "Wrote jets\n";
333    
334     return;
335     }
336    
337     //______________________________________________________________________________
338    
339     void MyVirtualJetProducer::inputTowers( )
340     {
341     std::vector<edm::Ptr<reco::Candidate> >::const_iterator inBegin = inputs_.begin(),
342     inEnd = inputs_.end(), i = inBegin;
343     for (; i != inEnd; ++i ) {
344     reco::CandidatePtr input = *i;
345     if (isnan(input->pt())) continue;
346     if (input->et() <inputEtMin_) continue;
347     if (input->energy()<inputEMin_) continue;
348     if (isAnomalousTower(input)) continue;
349 yilmaz 1.2
350     //Check consistency of kinematics
351     const CaloTower* ctc = dynamic_cast<const CaloTower*>(input.get());
352     if(ctc){
353     int ieta = ctc->id().ieta();
354     int iphi = ctc->id().iphi();
355    
356 yilmaz 1.6 if(0 && ntuple)ntuple->Fill(ieta, input->eta(), iphi, input->phi(),input->et(),ctc->emEt(),ctc->hadEt());
357 yilmaz 1.2
358     if(abs(ieta) < 5){
359    
360     if(0){
361     math::RhoEtaPhiVector v(1.4,input->eta(),input->phi());
362     GlobalPoint point(v.x(),v.y(),v.z());
363     // const DetId d = geo->getClosestCell(point);
364     // HcalDetId hd(d);
365     HcalDetId hd(0);
366     if(hd.ieta() != ieta || hd.iphi() != iphi){
367     cout<<"Inconsistent kinematics!!! ET = "<<input->pt()<<endl;
368     cout<<"ieta candidate : "<<ieta<<" ieta detid : "<<hd.ieta()<<endl;
369     cout<<"iphi candidate : "<<iphi<<" iphi detid : "<<hd.iphi()<<endl;
370     }
371    
372     }
373    
374     if(0){
375     HcalDetId det(HcalBarrel,ieta,iphi,1);
376    
377     if(geo->present(det)){
378     double eta = geo->getPosition(det).eta();
379     double phi = geo->getPosition(det).phi();
380    
381     if(input->eta() != eta || input->phi() != phi){
382     cout<<"Inconsistent kinematics!!! ET = "<<input->pt()<<endl;
383     cout<<"eta candidate : "<<input->eta()<<" eta detid : "<<eta<<endl;
384     cout<<"phi candidate : "<<input->phi()<<" phi detid : "<<phi<<endl;
385     }
386     }else{
387     cout<<"DetId not present in the Calo Geometry : ieta = "<<ieta<<" iphi = "<<iphi<<endl;
388     }
389     }
390     }
391    
392     }
393    
394 yilmaz 1.1 if (makeCaloJet(jetTypeE)&&doPVCorrection_) {
395     const CaloTower* tower=dynamic_cast<const CaloTower*>(input.get());
396     math::PtEtaPhiMLorentzVector ct(tower->p4(vertex_));
397     fjInputs_.push_back(fastjet::PseudoJet(ct.px(),ct.py(),ct.pz(),ct.energy()));
398     }
399     else {
400     fjInputs_.push_back(fastjet::PseudoJet(input->px(),input->py(),input->pz(),
401     input->energy()));
402     }
403     fjInputs_.back().set_user_index(i - inBegin);
404     }
405    
406     if ( restrictInputs_ && fjInputs_.size() > maxInputs_ ) {
407     reco::helper::GreaterByPtPseudoJet pTComparator;
408     std::sort(fjInputs_.begin(), fjInputs_.end(), pTComparator);
409     fjInputs_.resize(maxInputs_);
410     edm::LogWarning("JetRecoTooManyEntries") << "Too many inputs in the event, limiting to first " << maxInputs_ << ". Output is suspect.";
411     }
412     }
413    
414     //______________________________________________________________________________
415     bool MyVirtualJetProducer::isAnomalousTower(reco::CandidatePtr input)
416     {
417     if (!makeCaloJet(jetTypeE)) return false;
418     const CaloTower* tower=dynamic_cast<const CaloTower*>(input.get());
419     if (0==tower) return false;
420     if (tower->numBadEcalCells() >maxBadEcalCells_ ||
421     tower->numRecoveredEcalCells() >maxRecoveredEcalCells_ ||
422     tower->numProblematicEcalCells()>maxProblematicEcalCells_||
423     tower->numBadHcalCells() >maxBadHcalCells_ ||
424     tower->numRecoveredHcalCells() >maxRecoveredHcalCells_ ||
425     tower->numProblematicHcalCells()>maxProblematicHcalCells_) return true;
426     return false;
427     }
428    
429    
430     //------------------------------------------------------------------------------
431     // This is pure virtual.
432     //______________________________________________________________________________
433     // void MyVirtualJetProducer::runAlgorithm( edm::Event & iEvent, edm::EventSetup const& iSetup,
434     // std::vector<edm::Ptr<reco::Candidate> > const & inputs_);
435    
436     //______________________________________________________________________________
437     void MyVirtualJetProducer::copyConstituents(const vector<fastjet::PseudoJet>& fjConstituents,
438     reco::Jet* jet)
439     {
440     for (unsigned int i=0;i<fjConstituents.size();++i)
441     jet->addDaughter(inputs_[fjConstituents[i].user_index()]);
442     }
443    
444    
445     //______________________________________________________________________________
446     vector<reco::CandidatePtr>
447     MyVirtualJetProducer::getConstituents(const vector<fastjet::PseudoJet>&fjConstituents)
448     {
449     vector<reco::CandidatePtr> result;
450     for (unsigned int i=0;i<fjConstituents.size();i++) {
451     int index = fjConstituents[i].user_index();
452     reco::CandidatePtr candidate = inputs_[index];
453     result.push_back(candidate);
454     }
455     return result;
456     }
457    
458     void MyVirtualJetProducer::output(edm::Event & iEvent, edm::EventSetup const& iSetup)
459     {
460     // Write jets and constitutents. Will use fjJets_, inputs_
461     // and fjClusterSeq_
462     switch( jetTypeE ) {
463     case JetType::CaloJet :
464     writeJets<reco::CaloJet>( iEvent, iSetup);
465     break;
466     case JetType::PFJet :
467     writeJets<reco::PFJet>( iEvent, iSetup);
468     break;
469     case JetType::GenJet :
470     writeJets<reco::GenJet>( iEvent, iSetup);
471     break;
472     case JetType::TrackJet :
473     writeJets<reco::TrackJet>( iEvent, iSetup);
474     break;
475     case JetType::BasicJet :
476     writeJets<reco::BasicJet>( iEvent, iSetup);
477     break;
478     default:
479     throw cms::Exception("InvalidInput") << "invalid jet type in MyVirtualJetProducer\n";
480     break;
481     };
482    
483     }
484    
485     template< typename T >
486     void MyVirtualJetProducer::writeJets( edm::Event & iEvent, edm::EventSetup const& iSetup )
487     {
488     // produce output jet collection
489    
490     using namespace reco;
491    
492     std::auto_ptr<std::vector<T> > jets(new std::vector<T>() );
493     jets->reserve(fjJets_.size());
494    
495     for (unsigned int ijet=0;ijet<fjJets_.size();++ijet) {
496     // allocate this jet
497     T jet;
498     // get the fastjet jet
499     const fastjet::PseudoJet& fjJet = fjJets_[ijet];
500     // get the constituents from fastjet
501     std::vector<fastjet::PseudoJet> fjConstituents =
502     sorted_by_pt(fjClusterSeq_->constituents(fjJet));
503     // convert them to CandidatePtr vector
504     std::vector<CandidatePtr> constituents =
505     getConstituents(fjConstituents);
506    
507     // calcuate the jet area
508     double jetArea=0.0;
509     if ( doAreaFastjet_ ) {
510     fastjet::ClusterSequenceArea const * clusterSequenceWithArea =
511     dynamic_cast<fastjet::ClusterSequenceArea const *>(&*fjClusterSeq_);
512     jetArea = clusterSequenceWithArea->area(fjJet);
513     }
514    
515     // write the specifics to the jet (simultaneously sets 4-vector, vertex).
516     // These are overridden functions that will call the appropriate
517     // specific allocator.
518     writeSpecific(jet,
519     Particle::LorentzVector(fjJet.px(),
520     fjJet.py(),
521     fjJet.pz(),
522     fjJet.E()),
523     vertex_,
524     constituents, iSetup);
525    
526     jet.setJetArea (jetArea);
527    
528     if(doPUOffsetCorr_){
529     jet.setPileup(subtractor_->getPileUpEnergy(ijet));
530     }else{
531     jet.setPileup (0.0);
532     }
533    
534     // add to the list
535     jets->push_back(jet);
536     }
537    
538     // put the jets in the collection
539     iEvent.put(jets);
540    
541     // calculate rho (median pT per unit area, for PU&UE subtraction down the line
542     if (doRhoFastjet_) {
543 yilmaz 1.4 if(doFastJetNonUniform_){
544     std::auto_ptr<std::vector<double> > rhos(new std::vector<double>);
545     std::auto_ptr<std::vector<double> > sigmas(new std::vector<double>);
546     int nEta = puCenters_.size();
547     rhos->reserve(nEta);
548     sigmas->reserve(nEta);
549     fastjet::ClusterSequenceArea const * clusterSequenceWithArea =
550     dynamic_cast<fastjet::ClusterSequenceArea const *> ( &*fjClusterSeq_ );
551     for(int ie = 0; ie < nEta; ++ie){
552     double eta = puCenters_[ie];
553     double rho = 0;
554     double sigma = 0;
555     double etamin=eta-puWidth_;
556     double etamax=eta+puWidth_;
557     fastjet::RangeDefinition range_rho(etamin,etamax);
558     clusterSequenceWithArea->get_median_rho_and_sigma(range_rho, true, rho, sigma);
559     rhos->push_back(rho);
560     sigmas->push_back(sigma);
561     }
562     iEvent.put(rhos,"rhos");
563     iEvent.put(sigmas,"sigmas");
564     }else{
565     std::auto_ptr<double> rho(new double(0.0));
566     std::auto_ptr<double> sigma(new double(0.0));
567     double mean_area = 0;
568     fastjet::ClusterSequenceArea const * clusterSequenceWithArea =
569     dynamic_cast<fastjet::ClusterSequenceArea const *> ( &*fjClusterSeq_ );
570     clusterSequenceWithArea->get_median_rho_and_sigma(*fjRangeDef_,false,*rho,*sigma,mean_area);
571     iEvent.put(rho,"rho");
572     iEvent.put(sigma,"sigma");
573     }
574 yilmaz 1.1 }
575     }
576    
577 yilmaz 1.4