19 reRunAlgo_ (iConfig.getUntrackedParameter<bool>(
"reRunAlgo",
false)),
20 doAreaFastjet_ (iConfig.getParameter<bool> (
"doAreaFastjet")),
21 doRhoFastjet_ (iConfig.getParameter<bool> (
"doRhoFastjet")),
22 jetPtMin_(iConfig.getParameter<double> (
"jetPtMin")),
23 nSigmaPU_(iConfig.getParameter<double>(
"nSigmaPU")),
24 radiusPU_(iConfig.getParameter<double>(
"radiusPU")),
29 double ghostEtaMax = iConfig.
getParameter<
double>(
"Ghost_EtaMax");
31 int activeAreaRepeats = iConfig.
getParameter<
int> (
"Active_Area_Repeats");
33 double ghostArea = iConfig.
getParameter<
double> (
"GhostArea");
42 std::vector<fastjet::PseudoJet>& towers,
43 std::vector<fastjet::PseudoJet>&
output){
62 LogDebug(
"PileUpSubtractor")<<
"The subtractor setting up geometry...\n";
73 for(std::vector<DetId>::const_iterator did=alldid.begin(); did != alldid.end(); did++){
76 if( (hid).depth() == 1 ) {
79 if((hid).
ieta() != ietaold){
80 ietaold = (hid).
ieta();
103 LogDebug(
"PileUpSubtractor")<<
"The subtractor calculating pedestals...\n";
105 map<int,double> emean2;
106 map<int,int> ntowers;
108 int ietaold = -10000;
121 for (vector<fastjet::PseudoJet>::const_iterator input_object = coll.begin (),
122 fjInputsEnd = coll.end();
123 input_object != fjInputsEnd; ++input_object) {
126 ieta0 =
ieta( originalTower );
127 double Original_Et = originalTower->et();
129 if( ieta0-ietaold != 0 )
132 emean2[ieta0] = emean2[ieta0]+Original_Et*Original_Et;
139 emean2[ieta0] = emean2[ieta0]+Original_Et*Original_Et;
148 int it = (*gt).first;
150 double e1 = (*(
emean_.find(it))).second;
151 double e2 = (*emean2.find(it)).
second;
154 LogDebug(
"PileUpSubtractor")<<
" ieta : "<<it<<
" number of towers : "<<nt<<
" e1 : "<<e1<<
" e2 : "<<e2<<
"\n";
158 double eee = e2/nt - e1*e1/(nt*
nt);
178 LogDebug(
"PileUpSubtractor")<<
"The subtractor subtracting pedestals...\n";
183 for (vector<fastjet::PseudoJet>::iterator input_object = coll.begin (),
184 fjInputsEnd = coll.end();
185 input_object != fjInputsEnd; ++input_object) {
192 double etnew = itow->et() - (*(
emean_.find(it))).second - (*(
esigma_.find(it))).second;
193 float mScale = etnew/input_object->Et();
194 if(etnew < 0.) mScale = 0.;
197 input_object->pz()*mScale, input_object->e()*mScale);
199 int index = input_object->user_index();
200 input_object->reset ( towP4.px(),
204 input_object->set_user_index(index);
211 LogDebug(
"PileUpSubtractor")<<
"The subtractor calculating orphan input...\n";
215 vector<int> jettowers;
216 vector<pair<int,int> > excludedTowers;
218 vector <fastjet::PseudoJet>::iterator pseudojetTMP =
fjJets_->begin (),
221 for (; pseudojetTMP != fjJetsEnd ; ++pseudojetTMP) {
223 vector<fastjet::PseudoJet> newtowers;
230 excludedTowers.push_back(pair<int,int>(im->ieta(),im->iphi()));
234 vector<fastjet::PseudoJet>::const_iterator it =
fjInputs_->begin(),
237 for (; it != fjInputsEnd; ++it ) {
239 int index = it->user_index();
243 vector<pair<int,int> >::const_iterator exclude =
find(excludedTowers.begin(),excludedTowers.end(),pair<int,int>(ie,ip));
244 if(exclude != excludedTowers.end()) {
245 newtowers.push_back(*it);
246 jettowers.push_back(index);
254 for(vector<fastjet::PseudoJet>::const_iterator it =
fjInputs_->begin(),
255 fjInputsEnd =
fjInputs_->end(); it != fjInputsEnd; ++it ) {
257 int index = it->user_index();
258 vector<int>::const_iterator itjet =
find(jettowers.begin(),jettowers.end(),
index);
259 if( itjet == jettowers.end() ){
262 fastjet::PseudoJet orphan(originalTower->px(),originalTower->py(),originalTower->pz(),originalTower->energy());
263 orphan.set_user_index(index);
265 orphanInput.push_back(orphan);
274 LogDebug(
"PileUpSubtractor")<<
"The subtractor correcting jets...\n";
276 using namespace reco;
284 vector<fastjet::PseudoJet>::iterator pseudojetTMP =
fjJets_->begin (),
286 for (; pseudojetTMP != jetsEnd; ++pseudojetTMP) {
288 int ijet = pseudojetTMP -
fjJets_->begin();
291 std::vector<fastjet::PseudoJet> towers =
294 double newjetet = 0.;
295 for(vector<fastjet::PseudoJet>::const_iterator ito = towers.begin(),
296 towEnd = towers.end();
301 int it =
ieta( originalTower );
302 double Original_Et = originalTower->et();
304 if(etnew < 0.) etnew = 0;
305 newjetet = newjetet + etnew;
309 double mScale = newjetet/pseudojetTMP->Et();
310 LogDebug(
"PileUpSubtractor")<<
"pseudojetTMP->Et() : "<<pseudojetTMP->Et()<<
"\n";
311 LogDebug(
"PileUpSubtractor")<<
"newjetet : "<<newjetet<<
"\n";
313 LogDebug(
"PileUpSubtractor")<<
"pseudojetTMP->Et() - jetOffset_[ijet] : "<<pseudojetTMP->Et() -
jetOffset_[ijet]<<
"\n";
314 LogDebug(
"PileUpSubtractor")<<
"Scale is : "<<mScale<<
"\n";
316 int cshist = pseudojetTMP->cluster_hist_index();
317 pseudojetTMP->reset(pseudojetTMP->px()*mScale, pseudojetTMP->py()*mScale,
318 pseudojetTMP->pz()*mScale, pseudojetTMP->e()*mScale);
319 pseudojetTMP->set_cluster_hist_index(cshist);
346 throw cms::Exception(
"Invalid Constituent") <<
"CaloJet constituent is not of CaloTower type";
358 throw cms::Exception(
"Invalid Constituent") <<
"CaloJet constituent is not of CaloTower type";
T getParameter(std::string const &) const
virtual void calculateOrphanInput(std::vector< fastjet::PseudoJet > &orphanInput)
std::vector< double > jetOffset_
virtual double getMeanAtTower(const reco::CandidatePtr &in) const
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
boost::shared_ptr< fastjet::ClusterSequence > ClusterSequencePtr
std::vector< fastjet::PseudoJet > * fjJets_
std::map< int, double > esigma_
std::map< int, int > geomtowers_
std::vector< fastjet::PseudoJet > fjOriginalInputs_
virtual double getPileUpAtTower(const reco::CandidatePtr &in) const
PileUpSubtractor(const edm::ParameterSet &iConfig)
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
double deltaR(double eta1, double phi1, double eta2, double phi2)
int ieta(const reco::CandidatePtr &in) const
U second(std::pair< T, U > const &p)
virtual void offsetCorrectJets()
std::vector< fastjet::PseudoJet > * fjInputs_
std::map< int, int > ntowersWithJets_
virtual void setAlgorithm(ClusterSequencePtr &algorithm)
T const * get() const
Returns C++ pointer to the item.
int iphi() const
get the tower iphi
const GlobalPoint & getPosition(const DetId &id) const
Get the position of a given detector id.
virtual void subtractPedestal(std::vector< fastjet::PseudoJet > &coll)
#define EDM_REGISTER_PLUGINFACTORY(_factory_, _category_)
int iphi(const reco::CandidatePtr &in) const
CaloTowerDetId id() const
ActiveAreaSpecPtr fjActiveArea_
virtual double getSigmaAtTower(const reco::CandidatePtr &in) const
ClusterSequencePtr fjClusterSeq_
boost::shared_ptr< fastjet::RangeDefinition > RangeDefPtr
T const * product() const
boost::shared_ptr< fastjet::ActiveAreaSpec > ActiveAreaSpecPtr
std::vector< DetId > getValidDetIds() const
Get the list of all valid detector ids.
CaloGeometry const * geo_
std::map< int, double > emean_
virtual void calculatePedestal(std::vector< fastjet::PseudoJet > const &coll)
virtual void setupGeometryMap(edm::Event &iEvent, const edm::EventSetup &iSetup)
int ieta() const
get the tower ieta
virtual void reset(std::vector< edm::Ptr< reco::Candidate > > &input, std::vector< fastjet::PseudoJet > &towers, std::vector< fastjet::PseudoJet > &output)
std::vector< HcalDetId > allgeomid_
std::vector< edm::Ptr< reco::Candidate > > * inputs_