CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PileUpSubtractor.cc
Go to the documentation of this file.
1 
3 
8 
14 
15 #include <map>
16 using namespace std;
17 
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")),
25  geo_(0)
26 {
27  if ( doAreaFastjet_ || doRhoFastjet_ ) {
28  // default Ghost_EtaMax should be 5
29  double ghostEtaMax = iConfig.getParameter<double>("Ghost_EtaMax");
30  // default Active_Area_Repeats 1
31  int activeAreaRepeats = iConfig.getParameter<int> ("Active_Area_Repeats");
32  // default GhostArea 0.01
33  double ghostArea = iConfig.getParameter<double> ("GhostArea");
34  fjActiveArea_ = ActiveAreaSpecPtr(new fastjet::ActiveAreaSpec(ghostEtaMax,
35  activeAreaRepeats,
36  ghostArea));
37  fjRangeDef_ = RangeDefPtr( new fastjet::RangeDefinition(ghostEtaMax) );
38  }
39 }
40 
42  std::vector<fastjet::PseudoJet>& towers,
43  std::vector<fastjet::PseudoJet>& output){
44 
45  inputs_ = &input;
46  fjInputs_ = &towers;
47  fjJets_ = &output;
48  fjOriginalInputs_ = (*fjInputs_);
49  for (unsigned int i = 0; i < fjInputs_->size(); ++i){
50  fjOriginalInputs_[i].set_user_index((*fjInputs_)[i].user_index());
51  }
52 
53 }
54 
57 }
58 
60 {
61 
62  LogDebug("PileUpSubtractor")<<"The subtractor setting up geometry...\n";
63 
64  if(geo_ == 0) {
66  iSetup.get<CaloGeometryRecord>().get(pG);
67  geo_ = pG.product();
68  std::vector<DetId> alldid = geo_->getValidDetIds();
69 
70  int ietaold = -10000;
71  ietamax_ = -10000;
72  ietamin_ = 10000;
73  for(std::vector<DetId>::const_iterator did=alldid.begin(); did != alldid.end(); did++){
74  if( (*did).det() == DetId::Hcal ){
75  HcalDetId hid = HcalDetId(*did);
76  if( (hid).depth() == 1 ) {
77  allgeomid_.push_back(*did);
78 
79  if((hid).ieta() != ietaold){
80  ietaold = (hid).ieta();
81  geomtowers_[(hid).ieta()] = 1;
82  if((hid).ieta() > ietamax_) ietamax_ = (hid).ieta();
83  if((hid).ieta() < ietamin_) ietamin_ = (hid).ieta();
84  }
85  else{
86  geomtowers_[(hid).ieta()]++;
87  }
88  }
89  }
90  }
91  }
92 
93  for (int i = ietamin_; i<ietamax_+1; i++) {
94  ntowersWithJets_[i] = 0;
95  }
96 }
97 
98 //
99 // Calculate mean E and sigma from jet collection "coll".
100 //
101 void PileUpSubtractor::calculatePedestal( vector<fastjet::PseudoJet> const & coll )
102 {
103  LogDebug("PileUpSubtractor")<<"The subtractor calculating pedestals...\n";
104 
105  map<int,double> emean2;
106  map<int,int> ntowers;
107 
108  int ietaold = -10000;
109  int ieta0 = -100;
110 
111  // Initial values for emean_, emean2, esigma_, ntowers
112 
113  for(int i = ietamin_; i < ietamax_+1; i++)
114  {
115  emean_[i] = 0.;
116  emean2[i] = 0.;
117  esigma_[i] = 0.;
118  ntowers[i] = 0;
119  }
120 
121  for (vector<fastjet::PseudoJet>::const_iterator input_object = coll.begin (),
122  fjInputsEnd = coll.end();
123  input_object != fjInputsEnd; ++input_object) {
124 
125  const reco::CandidatePtr & originalTower=(*inputs_)[ input_object->user_index()];
126  ieta0 = ieta( originalTower );
127  double Original_Et = originalTower->et();
128 
129  if( ieta0-ietaold != 0 )
130  {
131  emean_[ieta0] = emean_[ieta0]+Original_Et;
132  emean2[ieta0] = emean2[ieta0]+Original_Et*Original_Et;
133  ntowers[ieta0] = 1;
134  ietaold = ieta0;
135  }
136  else
137  {
138  emean_[ieta0] = emean_[ieta0]+Original_Et;
139  emean2[ieta0] = emean2[ieta0]+Original_Et*Original_Et;
140  ntowers[ieta0]++;
141  }
142 
143  }
144 
145  for(map<int,int>::const_iterator gt = geomtowers_.begin(); gt != geomtowers_.end(); gt++)
146  {
147 
148  int it = (*gt).first;
149 
150  double e1 = (*(emean_.find(it))).second;
151  double e2 = (*emean2.find(it)).second;
152  int nt = (*gt).second - (*(ntowersWithJets_.find(it))).second;
153 
154  LogDebug("PileUpSubtractor")<<" ieta : "<<it<<" number of towers : "<<nt<<" e1 : "<<e1<<" e2 : "<<e2<<"\n";
155 
156  if(nt > 0) {
157  emean_[it] = e1/nt;
158  double eee = e2/nt - e1*e1/(nt*nt);
159  if(eee<0.) eee = 0.;
160  esigma_[it] = nSigmaPU_*sqrt(eee);
161  }
162  else
163  {
164  emean_[it] = 0.;
165  esigma_[it] = 0.;
166  }
167  LogDebug("PileUpSubtractor")<<" ieta : "<<it<<" Pedestals : "<<emean_[it]<<" "<<esigma_[it]<<"\n";
168  }
169 }
170 
171 
172 //
173 // Subtract mean and sigma from fjInputs_
174 //
175 void PileUpSubtractor::subtractPedestal(vector<fastjet::PseudoJet> & coll)
176 {
177 
178  LogDebug("PileUpSubtractor")<<"The subtractor subtracting pedestals...\n";
179 
180  int it = -100;
181  int ip = -100;
182 
183  for (vector<fastjet::PseudoJet>::iterator input_object = coll.begin (),
184  fjInputsEnd = coll.end();
185  input_object != fjInputsEnd; ++input_object) {
186 
187  reco::CandidatePtr const & itow = (*inputs_)[ input_object->user_index() ];
188 
189  it = ieta( itow );
190  ip = iphi( itow );
191 
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.;
195 
196  math::XYZTLorentzVectorD towP4(input_object->px()*mScale, input_object->py()*mScale,
197  input_object->pz()*mScale, input_object->e()*mScale);
198 
199  int index = input_object->user_index();
200  input_object->reset ( towP4.px(),
201  towP4.py(),
202  towP4.pz(),
203  towP4.energy() );
204  input_object->set_user_index(index);
205  }
206 }
207 
208 void PileUpSubtractor::calculateOrphanInput(vector<fastjet::PseudoJet> & orphanInput)
209 {
210 
211  LogDebug("PileUpSubtractor")<<"The subtractor calculating orphan input...\n";
212 
213  (*fjInputs_) = fjOriginalInputs_;
214 
215  vector<int> jettowers; // vector of towers indexed by "user_index"
216  vector<pair<int,int> > excludedTowers; // vector of excluded ieta, iphi values
217 
218  vector <fastjet::PseudoJet>::iterator pseudojetTMP = fjJets_->begin (),
219  fjJetsEnd = fjJets_->end();
220 
221  for (; pseudojetTMP != fjJetsEnd ; ++pseudojetTMP) {
222 
223  vector<fastjet::PseudoJet> newtowers;
224  // find towers within radiusPU_ of this jet
225  for(vector<HcalDetId>::const_iterator im = allgeomid_.begin(); im != allgeomid_.end(); im++)
226  {
227  double dr = reco::deltaR(geo_->getPosition((DetId)(*im)),(*pseudojetTMP));
228  if( dr < radiusPU_) {
229  ntowersWithJets_[(*im).ieta()]++;
230  excludedTowers.push_back(pair<int,int>(im->ieta(),im->iphi()));
231  }
232  }
233 
234  vector<fastjet::PseudoJet>::const_iterator it = fjInputs_->begin(),
235  fjInputsEnd = fjInputs_->end();
236 
237  for (; it != fjInputsEnd; ++it ) {
238 
239  int index = it->user_index();
240  int ie = ieta((*inputs_)[index]);
241  int ip = iphi((*inputs_)[index]);
242 
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);
247  } //dr < 0.5
248  } // initial input collection
249  } // pseudojets
250 
251  //
252  // Create a new collections from the towers not included in jets
253  //
254  for(vector<fastjet::PseudoJet>::const_iterator it = fjInputs_->begin(),
255  fjInputsEnd = fjInputs_->end(); it != fjInputsEnd; ++it ) {
256 
257  int index = it->user_index();
258  vector<int>::const_iterator itjet = find(jettowers.begin(),jettowers.end(),index);
259  if( itjet == jettowers.end() ){
260 
261  const reco::CandidatePtr& originalTower = (*inputs_)[index];
262  fastjet::PseudoJet orphan(originalTower->px(),originalTower->py(),originalTower->pz(),originalTower->energy());
263  orphan.set_user_index(index);
264 
265  orphanInput.push_back(orphan);
266  }
267  }
268 }
269 
270 
272 {
273 
274  LogDebug("PileUpSubtractor")<<"The subtractor correcting jets...\n";
275  jetOffset_.clear();
276  using namespace reco;
277 
278  //
279  // Reestimate energy of jet (energy of jet with initial map)
280  //
281 
282  jetOffset_.reserve(fjJets_->size());
283 
284  vector<fastjet::PseudoJet>::iterator pseudojetTMP = fjJets_->begin (),
285  jetsEnd = fjJets_->end();
286  for (; pseudojetTMP != jetsEnd; ++pseudojetTMP) {
287 
288  int ijet = pseudojetTMP - fjJets_->begin();
289  jetOffset_[ijet] = 0;
290 
291  std::vector<fastjet::PseudoJet> towers =
292  sorted_by_pt(fjClusterSeq_->constituents(*pseudojetTMP));
293 
294  double newjetet = 0.;
295  for(vector<fastjet::PseudoJet>::const_iterator ito = towers.begin(),
296  towEnd = towers.end();
297  ito != towEnd;
298  ++ito)
299  {
300  const reco::CandidatePtr& originalTower = (*inputs_)[ito->user_index()];
301  int it = ieta( originalTower );
302  double Original_Et = originalTower->et();
303  double etnew = Original_Et - (*emean_.find(it)).second - (*esigma_.find(it)).second;
304  if(etnew < 0.) etnew = 0;
305  newjetet = newjetet + etnew;
306  jetOffset_[ijet] += Original_Et - etnew;
307  }
308 
309  double mScale = newjetet/pseudojetTMP->Et();
310  LogDebug("PileUpSubtractor")<<"pseudojetTMP->Et() : "<<pseudojetTMP->Et()<<"\n";
311  LogDebug("PileUpSubtractor")<<"newjetet : "<<newjetet<<"\n";
312  LogDebug("PileUpSubtractor")<<"jetOffset_[ijet] : "<<jetOffset_[ijet]<<"\n";
313  LogDebug("PileUpSubtractor")<<"pseudojetTMP->Et() - jetOffset_[ijet] : "<<pseudojetTMP->Et() - jetOffset_[ijet]<<"\n";
314  LogDebug("PileUpSubtractor")<<"Scale is : "<<mScale<<"\n";
315 
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);
320 
321  }
322 }
323 
325  int it = ieta(in);
326  return (*emean_.find(it)).second;
327 }
328 
330  int it = ieta(in);
331  return (*esigma_.find(it)).second;
332 }
333 
335  int it = ieta(in);
336  return (*emean_.find(it)).second + (*esigma_.find(it)).second;
337 }
338 
340  int it = 0;
341  const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
342  if(ctc){
343  it = ctc->id().ieta();
344  } else
345  {
346  throw cms::Exception("Invalid Constituent") << "CaloJet constituent is not of CaloTower type";
347  }
348  return it;
349 }
350 
352  int it = 0;
353  const CaloTower* ctc = dynamic_cast<const CaloTower*>(in.get());
354  if(ctc){
355  it = ctc->id().iphi();
356  } else
357  {
358  throw cms::Exception("Invalid Constituent") << "CaloJet constituent is not of CaloTower type";
359  }
360  return it;
361 }
362 
363 
365 EDM_REGISTER_PLUGINFACTORY(PileUpSubtractorFactory,"PileUpSubtractorFactory");
366 
367 
368 
#define LogDebug(id)
T getParameter(std::string const &) const
virtual void calculateOrphanInput(std::vector< fastjet::PseudoJet > &orphanInput)
int i
Definition: DBlmapReader.cc:9
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.
Definition: LorentzVector.h:15
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)
Definition: FindCaloHit.cc:7
double deltaR(double eta1, double phi1, double eta2, double phi2)
Definition: deltaR.h:19
int ieta(const reco::CandidatePtr &in) const
U second(std::pair< T, U > const &p)
virtual void offsetCorrectJets()
RangeDefPtr fjRangeDef_
std::vector< fastjet::PseudoJet > * fjInputs_
std::map< int, int > ntowersWithJets_
int iEvent
Definition: GenABIO.cc:243
virtual void setAlgorithm(ClusterSequencePtr &algorithm)
T const * get() const
Returns C++ pointer to the item.
Definition: Ptr.h:147
T sqrt(T t)
Definition: SSEVec.h:28
int iphi() const
get the tower iphi
const GlobalPoint & getPosition(const DetId &id) const
Get the position of a given detector id.
Definition: CaloGeometry.cc:68
virtual void subtractPedestal(std::vector< fastjet::PseudoJet > &coll)
int nt
Definition: AMPTWrapper.h:32
tuple input
Definition: collect_tpl.py:10
#define EDM_REGISTER_PLUGINFACTORY(_factory_, _category_)
Definition: DetId.h:20
int iphi(const reco::CandidatePtr &in) const
CaloTowerDetId id() const
Definition: CaloTower.h:72
JetCorrectorParametersCollection coll
Definition: classes.h:14
ActiveAreaSpecPtr fjActiveArea_
virtual double getSigmaAtTower(const reco::CandidatePtr &in) const
ClusterSequencePtr fjClusterSeq_
const T & get() const
Definition: EventSetup.h:55
boost::shared_ptr< fastjet::RangeDefinition > RangeDefPtr
T const * product() const
Definition: ESHandle.h:62
boost::shared_ptr< fastjet::ActiveAreaSpec > ActiveAreaSpecPtr
std::vector< DetId > getValidDetIds() const
Get the list of all valid detector ids.
Definition: CaloGeometry.cc:90
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)
DTDigi & gt
std::vector< HcalDetId > allgeomid_
std::vector< edm::Ptr< reco::Candidate > > * inputs_