CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DistortedPFCandProducer.cc
Go to the documentation of this file.
1 #include <memory>
5 
8 
9 //
10 // class declaration
11 //
13  public:
16 
17  private:
18  virtual void beginJob() ;
19  virtual void produce(edm::Event&, const edm::EventSetup&);
20  virtual void endJob() ;
21 
25  std::vector<double> etaBinEdges_;
26 
27  std::vector<double> shiftOnOneOverPt_; // in [1/GeV]
28  std::vector<double> relativeShiftOnPt_;
29  std::vector<double> uncertaintyOnOneOverPt_; // in [1/GeV]
30  std::vector<double> relativeUncertaintyOnPt_;
31 
32  std::vector<double> efficiencyRatioOverMC_;
33 };
34 
45 
46 #include <CLHEP/Random/RandFlat.h>
47 #include <CLHEP/Random/RandGauss.h>
48 
51 
52  // What is being produced
53  produces<std::vector<reco::PFCandidate> >();
54 
55  // Input products
56  muonTag_ = pset.getUntrackedParameter<edm::InputTag> ("MuonTag", edm::InputTag("muons"));
57  pfTag_ = pset.getUntrackedParameter<edm::InputTag> ("PFTag", edm::InputTag("particleFlow"));
58  genMatchMapTag_ = pset.getUntrackedParameter<edm::InputTag> ("GenMatchMapTag", edm::InputTag("genMatchMap"));
59 
60  // Eta edges
61  std::vector<double> defEtaEdges;
62  defEtaEdges.push_back(-999999.);
63  defEtaEdges.push_back(999999.);
64  etaBinEdges_ = pset.getUntrackedParameter<std::vector<double> > ("EtaBinEdges",defEtaEdges);
65  unsigned int ninputs_expected = etaBinEdges_.size()-1;
66 
67  // Distortions in muon momentum
68  std::vector<double> defDistortion;
69  defDistortion.push_back(0.);
70 
71  shiftOnOneOverPt_ = pset.getUntrackedParameter<std::vector<double> > ("ShiftOnOneOverPt",defDistortion); // in [1/GeV]
72  if (shiftOnOneOverPt_.size()==1 && ninputs_expected>1) {
73  for (unsigned int i=1; i<ninputs_expected; i++){ shiftOnOneOverPt_.push_back(shiftOnOneOverPt_[0]);}
74  }
75 
76  relativeShiftOnPt_ = pset.getUntrackedParameter<std::vector<double> > ("RelativeShiftOnPt",defDistortion);
77  if (relativeShiftOnPt_.size()==1 && ninputs_expected>1) {
78  for (unsigned int i=1; i<ninputs_expected; i++){ relativeShiftOnPt_.push_back(relativeShiftOnPt_[0]);}
79  }
80 
81  uncertaintyOnOneOverPt_ = pset.getUntrackedParameter<std::vector<double> > ("UncertaintyOnOneOverPt",defDistortion); // in [1/GeV]
82  if (uncertaintyOnOneOverPt_.size()==1 && ninputs_expected>1) {
83  for (unsigned int i=1; i<ninputs_expected; i++){ uncertaintyOnOneOverPt_.push_back(uncertaintyOnOneOverPt_[0]);}
84  }
85 
86  relativeUncertaintyOnPt_ = pset.getUntrackedParameter<std::vector<double> > ("RelativeUncertaintyOnPt",defDistortion);
87  if (relativeUncertaintyOnPt_.size()==1 && ninputs_expected>1) {
88  for (unsigned int i=1; i<ninputs_expected; i++){ relativeUncertaintyOnPt_.push_back(relativeUncertaintyOnPt_[0]);}
89  }
90 
91  // Data/MC efficiency ratios
92  std::vector<double> defEfficiencyRatio;
93  defEfficiencyRatio.push_back(1.);
94  efficiencyRatioOverMC_ = pset.getUntrackedParameter<std::vector<double> > ("EfficiencyRatioOverMC",defEfficiencyRatio);
95  if (efficiencyRatioOverMC_.size()==1 && ninputs_expected>1) {
96  for (unsigned int i=1; i<ninputs_expected; i++){ efficiencyRatioOverMC_.push_back(efficiencyRatioOverMC_[0]);}
97  }
98 
99  // Send a warning if there are inconsistencies in vector sizes !!
100  bool effWrong = efficiencyRatioOverMC_.size()!=ninputs_expected;
101  bool momWrong = shiftOnOneOverPt_.size()!=ninputs_expected
102  || relativeShiftOnPt_.size()!=ninputs_expected
103  || uncertaintyOnOneOverPt_.size()!=ninputs_expected
104  || relativeUncertaintyOnPt_.size()!=ninputs_expected;
105  if ( effWrong and momWrong) {
106  edm::LogError("") << "WARNING: DistortedPFCandProducer : Size of some parameters do not match the EtaBinEdges vector!!";
107  }
108 
109 }
110 
113 }
114 
117 }
118 
121 }
122 
125 
126  if (ev.isRealData()) return;
127 
128  // Muon collection
129  edm::Handle<edm::View<reco::Muon> > muonCollection;
130  if (!ev.getByLabel(muonTag_, muonCollection)) {
131  edm::LogError("") << ">>> Muon collection does not exist !!!";
132  return;
133  }
134 
136  if (!ev.getByLabel(genMatchMapTag_, genMatchMap)) {
137  edm::LogError("") << ">>> Muon-GenParticle match map does not exist !!!";
138  return;
139  }
140 
141 
142  // Get PFCandidate collection
144  if (!ev.getByLabel(pfTag_, pfCollection)) {
145  edm::LogError("") << ">>> PFCandidate collection does not exist !!!";
146  return;
147  }
148 
149  unsigned int muonCollectionSize = muonCollection->size();
150  unsigned int pfCollectionSize = pfCollection->size();
151 
152  if (pfCollectionSize<1) return;
153 
154 
155  // Ask for PfMuon consistency
156  bool pfMuonFound = false;
157 
158 
159 
160  std::auto_ptr<reco::PFCandidateCollection> newmuons (new reco::PFCandidateCollection);
161 
162 
163  // Loop on all PF candidates
164  for (unsigned int j=0; j<pfCollectionSize; j++) {
165  edm::RefToBase<reco::PFCandidate> pf = pfCollection->refAt(j);
166 
167 
168  // New PF muon
169  double ptmu = pf->pt();
170 
171 
172  for (unsigned int i=0; i<muonCollectionSize; i++) {
173  edm::RefToBase<reco::Muon> mu = muonCollection->refAt(i);
174 
175 
176  // Check the muon is in the PF collection
177  if (pf->particleId()==reco::PFCandidate::mu) {
178  reco::MuonRef muref = pf->muonRef();
179  if (muref.isNonnull()) {
180  if (muref.key()==mu.key()) {
181  if ( mu->isStandAloneMuon() && ptmu == muref->standAloneMuon()->pt() && (
182  ( !mu->isGlobalMuon() || ( mu->isGlobalMuon() && ptmu != muref->combinedMuon()->pt() ) ) &&
183  ( !mu->isTrackerMuon() || ( mu->isTrackerMuon() && ptmu != mu->track()->pt() ) ) )
184  ) {
185  pfMuonFound = false;
186  }
187  else if ( !mu->isTrackerMuon() ){
188  pfMuonFound = false;
189  }
190  else{
191  pfMuonFound = true;}
192  }
193  else {pfMuonFound = false; }
194 
195  }
196  }
197 
198  // do nothing if StandAlone muon
199  //const reco::Track& track = *trackRef;
200 
201  if ( !pfMuonFound) continue;
202 
203  double ptgen = pf->pt();
204  double etagen = pf->eta();
205 
206 
207  reco::GenParticleRef gen = (*genMatchMap)[mu];
208  if( !gen.isNull()) {
209  ptgen = gen->pt();
210  etagen = gen->eta();
211  LogTrace("") << ">>> Muon-GenParticle match found; ptmu= " << pf->pt() << ", ptgen= " << ptgen;
212  } else {
213  LogTrace("") << ">>> MUON-GENPARTICLE MATCH NOT FOUND!!!";
214  }
215 
216 
217  // Initialize parameters
218  double effRatio = 0.;
219  double shift1 = 0.;
220  double shift2 = 0.;
221  double sigma1 = 0.;
222  double sigma2 = 0.;
223 
224  // Find out which eta bin should be used
225  unsigned int nbins = etaBinEdges_.size()-1;
226  unsigned int etaBin = nbins;
227  if (etagen>etaBinEdges_[0] && etagen<etaBinEdges_[nbins]) {
228  for (unsigned int j=1; j<=nbins; ++j) {
229  if (etagen>etaBinEdges_[j]) continue;
230  etaBin = j-1;
231  break;
232  }
233  }
234  if (etaBin<nbins) {
235  LogTrace("") << ">>> etaBin: " << etaBin << ", for etagen =" << etagen;
236  } else {
237  // Muon is rejected if outside the considered eta range
238  LogTrace("") << ">>> Muon outside eta range: reject it; etagen = " << etagen;
239  pfMuonFound = false;
240  continue;
241  }
242 
243  if (!pfMuonFound) continue;
244 
245  // Set shifts
246  shift1 = shiftOnOneOverPt_[etaBin];
247  shift2 = relativeShiftOnPt_[etaBin];
248  LogTrace("") << "\tshiftOnOneOverPt= " << shift1*100 << " [%]";
249  LogTrace("") << "\trelativeShiftOnPt= " << shift2*100 << " [%]";
250 
251  // Set resolutions
252  sigma1 = uncertaintyOnOneOverPt_[etaBin];
253  sigma2 = relativeUncertaintyOnPt_[etaBin];
254  LogTrace("") << "\tuncertaintyOnOneOverPt= " << sigma1 << " [1/GeV]";
255  LogTrace("") << "\trelativeUncertaintyOnPt= " << sigma2*100 << " [%]";
256 
257  // Set efficiency ratio
258  effRatio = efficiencyRatioOverMC_[etaBin];
259  LogTrace("") << "\tefficiencyRatioOverMC= " << effRatio;
260 
261  // Reject muons according to efficiency ratio
262  double rndf = CLHEP::RandFlat::shoot();
263  if (rndf>effRatio) continue;
264 
265  // Gaussian Random numbers for smearing
266  double rndg1 = CLHEP::RandGauss::shoot();
267  double rndg2 = CLHEP::RandGauss::shoot();
268 
269 
270  // change here the pt of the candidate, if it is a muon
271 
272  ptmu += ptgen * ( shift1*ptgen + shift2 + sigma1*rndg1*ptgen + sigma2*rndg2);
273  pfMuonFound = false ;
274 
275  }
276 
277  reco::PFCandidate* newmu = pf->clone();
278  newmu->setP4 (
280  ptmu, pf->eta(), pf->phi(), pf->mass()
281  )
282  );
283 
284  newmuons->push_back(*newmu);
285  }
286 
287 
288  ev.put(newmuons);
289 
290 
291 }
292 
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
std::vector< double > relativeUncertaintyOnPt_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
virtual PFCandidate * clone() const
return a clone
Definition: PFCandidate.cc:118
bool isTrackerMuon() const
Definition: Muon.h:212
virtual TrackRef track() const
reference to a Track
Definition: Muon.h:50
bool isGlobalMuon() const
Definition: Muon.h:211
virtual void setP4(const LorentzVector &p4)
set 4-momentum
bool isStandAloneMuon() const
Definition: Muon.h:213
bool isRealData() const
Definition: EventBase.h:60
virtual double eta() const
momentum pseudorapidity
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
virtual double mass() const
mass
bool isNull() const
Checks for null.
Definition: Ref.h:247
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:85
std::vector< double > uncertaintyOnOneOverPt_
int j
Definition: DBlmapReader.cc:9
virtual void produce(edm::Event &, const edm::EventSetup &)
const int mu
Definition: Constants.h:23
size_t key() const
Definition: RefToBase.h:228
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
#define LogTrace(id)
reco::MuonRef muonRef() const
Definition: PFCandidate.cc:355
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
virtual double pt() const
transverse momentum
std::vector< double > efficiencyRatioOverMC_
key_type key() const
Accessor for product key.
Definition: Ref.h:266
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:33
math::PtEtaPhiMLorentzVector PolarLorentzVector
Lorentz vector.
Definition: Particle.h:27
virtual ParticleType particleId() const
Definition: PFCandidate.h:324
std::vector< double > relativeShiftOnPt_
std::vector< double > shiftOnOneOverPt_
virtual double phi() const
momentum azimuthal angle
std::vector< double > etaBinEdges_
DistortedPFCandProducer(const edm::ParameterSet &)