CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ZToMuMuGammaAnalyzer.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 //
4 
7 
8 
19 using namespace std;
20 
21 
23 {
24 
25  fName_ = pset.getUntrackedParameter<string>("Name");
26  verbosity_ = pset.getUntrackedParameter<int>("Verbosity");
27  prescaleFactor_ = pset.getUntrackedParameter<int>("prescaleFactor",1);
28  standAlone_ = pset.getParameter<bool>("standAlone");
29  outputFileName_ = pset.getParameter<string>("OutputFileName");
30  isHeavyIon_ = pset.getUntrackedParameter<bool>("isHeavyIon",false);
31  triggerEvent_ = pset.getParameter<edm::InputTag>("triggerEvent");
32  useTriggerFiltering_ = pset.getParameter<bool>("useTriggerFiltering");
33  //
34  photonProducer_ = pset.getParameter<string>("phoProducer");
35  photonCollection_ = pset.getParameter<string>("photonCollection");
36 
37  barrelRecHitProducer_ = pset.getParameter<string>("barrelRecHitProducer");
38  barrelRecHitCollection_ = pset.getParameter<string>("barrelRecHitCollection");
39 
40  endcapRecHitProducer_ = pset.getParameter<string>("endcapRecHitProducer");
41  endcapRecHitCollection_ = pset.getParameter<string>("endcapRecHitCollection");
42 
43  muonProducer_ = pset.getParameter<string>("muonProducer");
44  muonCollection_ = pset.getParameter<string>("muonCollection");
45  // Muon selection
46  muonMinPt_ = pset.getParameter<double>("muonMinPt");
47  minPixStripHits_ = pset.getParameter<int>("minPixStripHits");
48  muonMaxChi2_ = pset.getParameter<double>("muonMaxChi2");
49  muonMaxDxy_ = pset.getParameter<double>("muonMaxDxy");
50  muonMatches_ = pset.getParameter<int>("muonMatches");
51  validPixHits_ = pset.getParameter<int>("validPixHits");
52  validMuonHits_ = pset.getParameter<int>("validMuonHits");
53  muonTrackIso_ = pset.getParameter<double>("muonTrackIso");
54  muonTightEta_ = pset.getParameter<double>("muonTightEta");
55  // Dimuon selection
56  minMumuInvMass_ = pset.getParameter<double>("minMumuInvMass");
57  maxMumuInvMass_ = pset.getParameter<double>("maxMumuInvMass");
58  // Photon selection
59  photonMinEt_ = pset.getParameter<double>("photonMinEt");
60  photonMaxEta_ = pset.getParameter<double>("photonMaxEta");
61  photonTrackIso_ = pset.getParameter<double>("photonTrackIso");
62  // mumuGamma selection
63  nearMuonDr_ = pset.getParameter<double>("nearMuonDr");
64  nearMuonHcalIso_ = pset.getParameter<double>("nearMuonHcalIso");
65  farMuonEcalIso_ = pset.getParameter<double>("farMuonEcalIso");
66  farMuonTrackIso_ = pset.getParameter<double>("farMuonTrackIso");
67  farMuonMinPt_ = pset.getParameter<double>("farMuonMinPt");
68  minMumuGammaInvMass_ = pset.getParameter<double>("minMumuGammaInvMass");
69  maxMumuGammaInvMass_ = pset.getParameter<double>("maxMumuGammaInvMass");
70  //
71  parameters_ = pset;
72 
73 }
74 
75 
76 
78 
79 
81 {
82 
83  nEvt_=0;
84  nEntry_=0;
85 
86  dbe_ = 0;
88 
89 
90 
91 // double eMin = parameters_.getParameter<double>("eMin");
92 // double eMax = parameters_.getParameter<double>("eMax");
93 // int eBin = parameters_.getParameter<int>("eBin");
94 
95  double etMin = parameters_.getParameter<double>("etMin");
96  double etMax = parameters_.getParameter<double>("etMax");
97  int etBin = parameters_.getParameter<int>("etBin");
98 
99 // double sumMin = parameters_.getParameter<double>("sumMin");
100 // double sumMax = parameters_.getParameter<double>("sumMax");
101 // int sumBin = parameters_.getParameter<int>("sumBin");
102 
103 // double etaMin = parameters_.getParameter<double>("etaMin");
104 // double etaMax = parameters_.getParameter<double>("etaMax");
105 // int etaBin = parameters_.getParameter<int>("etaBin");
106 
107 // double phiMin = parameters_.getParameter<double>("phiMin");
108 // double phiMax = parameters_.getParameter<double>("phiMax");
109 // int phiBin = parameters_.getParameter<int>("phiBin");
110 
111 // double r9Min = parameters_.getParameter<double>("r9Min");
112 // double r9Max = parameters_.getParameter<double>("r9Max");
113 // int r9Bin = parameters_.getParameter<int>("r9Bin");
114 
115 // double hOverEMin = parameters_.getParameter<double>("hOverEMin");
116 // double hOverEMax = parameters_.getParameter<double>("hOverEMax");
117 // int hOverEBin = parameters_.getParameter<int>("hOverEBin");
118 
119 // double xMin = parameters_.getParameter<double>("xMin");
120 // double xMax = parameters_.getParameter<double>("xMax");
121 // int xBin = parameters_.getParameter<int>("xBin");
122 
123 // double yMin = parameters_.getParameter<double>("yMin");
124 // double yMax = parameters_.getParameter<double>("yMax");
125 // int yBin = parameters_.getParameter<int>("yBin");
126 
127 // double numberMin = parameters_.getParameter<double>("numberMin");
128 // double numberMax = parameters_.getParameter<double>("numberMax");
129 // int numberBin = parameters_.getParameter<int>("numberBin");
130 
131 // double zMin = parameters_.getParameter<double>("zMin");
132 // double zMax = parameters_.getParameter<double>("zMax");
133 // int zBin = parameters_.getParameter<int>("zBin");
134 
135 // double rMin = parameters_.getParameter<double>("rMin");
136 // double rMax = parameters_.getParameter<double>("rMax");
137 // int rBin = parameters_.getParameter<int>("rBin");
138 
139 // double dPhiTracksMin = parameters_.getParameter<double>("dPhiTracksMin");
140 // double dPhiTracksMax = parameters_.getParameter<double>("dPhiTracksMax");
141 // int dPhiTracksBin = parameters_.getParameter<int>("dPhiTracksBin");
142 
143 // double dEtaTracksMin = parameters_.getParameter<double>("dEtaTracksMin");
144 // double dEtaTracksMax = parameters_.getParameter<double>("dEtaTracksMax");
145 // int dEtaTracksBin = parameters_.getParameter<int>("dEtaTracksBin");
146 
147 // double sigmaIetaMin = parameters_.getParameter<double>("sigmaIetaMin");
148 // double sigmaIetaMax = parameters_.getParameter<double>("sigmaIetaMax");
149 // int sigmaIetaBin = parameters_.getParameter<int>("sigmaIetaBin");
150 
151 // double eOverPMin = parameters_.getParameter<double>("eOverPMin");
152 // double eOverPMax = parameters_.getParameter<double>("eOverPMax");
153 // int eOverPBin = parameters_.getParameter<int>("eOverPBin");
154 
155 // double chi2Min = parameters_.getParameter<double>("chi2Min");
156 // double chi2Max = parameters_.getParameter<double>("chi2Max");
157 // int chi2Bin = parameters_.getParameter<int>("chi2Bin");
158 
159 
160 // int reducedEtBin = etBin/4;
161 // int reducedEtaBin = etaBin/4;
162 // int reducedSumBin = sumBin/4;
163 // int reducedR9Bin = r9Bin/4;
164 
165 
167 
168  if (dbe_) {
169 
170 
171  dbe_->setCurrentFolder("Egamma/PhotonAnalyzer/ZToMuMuGamma");
172  h1_mumuInvMass_ = dbe_->book1D("mumuInvMass","Two muon invariant mass: M (GeV)",etBin/2,etMin,etMax/2);
173  h1_mumuGammaInvMass_ = dbe_->book1D("mumuGammaInvMass","Two-muon plus gamma invariant mass: M (GeV)",etBin/2,etMin,etMax/2);
174 
175 
176 
177  }//end if(dbe_)
178 
179 
180 }//end BeginJob
181 
182 
183 
185 {
186  using namespace edm;
187 
188  if (nEvt_% prescaleFactor_ ) return;
189  nEvt_++;
190  LogInfo("ZToMuMuGammaAnalyzer") << "ZToMuMuGammaAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
191 
192 
193  // Get the trigger results
194  bool validTriggerEvent=true;
195  edm::Handle<trigger::TriggerEvent> triggerEventHandle;
196  trigger::TriggerEvent triggerEvent;
197  e.getByLabel(triggerEvent_,triggerEventHandle);
198  if(!triggerEventHandle.isValid()) {
199  edm::LogInfo("PhotonAnalyzer") << "Error! Can't get the product "<< triggerEvent_.label() << endl;
200  validTriggerEvent=false;
201  }
202  if(validTriggerEvent) triggerEvent = *(triggerEventHandle.product());
203 
204 
205  // Get the reconstructed photons
206  bool validPhotons=true;
207  Handle<reco::PhotonCollection> photonHandle;
208  reco::PhotonCollection photonCollection;
209  e.getByLabel(photonProducer_, photonCollection_ , photonHandle);
210  if ( !photonHandle.isValid()) {
211  edm::LogInfo("ZToMuMuGammaAnalyzer") << "Error! Can't get the product "<< photonCollection_ << endl;
212  validPhotons=false;
213  }
214  if(validPhotons) photonCollection = *(photonHandle.product());
215 
216  // Get the PhotonId objects
217  bool validloosePhotonID=true;
218  Handle<edm::ValueMap<bool> > loosePhotonFlag;
219  edm::ValueMap<bool> loosePhotonID;
220  e.getByLabel("PhotonIDProd", "PhotonCutBasedIDLoose", loosePhotonFlag);
221  if ( !loosePhotonFlag.isValid()) {
222  edm::LogInfo("ZToMuMuGammaAnalyzer") << "Error! Can't get the product "<< "PhotonCutBasedIDLoose" << endl;
223  validloosePhotonID=false;
224  }
225  if (validloosePhotonID) loosePhotonID = *(loosePhotonFlag.product());
226 
227  bool validtightPhotonID=true;
228  Handle<edm::ValueMap<bool> > tightPhotonFlag;
229  edm::ValueMap<bool> tightPhotonID;
230  e.getByLabel("PhotonIDProd", "PhotonCutBasedIDTight", tightPhotonFlag);
231  if ( !tightPhotonFlag.isValid()) {
232  edm::LogInfo("ZToMuMuGammaAnalyzer") << "Error! Can't get the product "<< "PhotonCutBasedIDTight" << endl;
233  validtightPhotonID=false;
234  }
235  if (validtightPhotonID) tightPhotonID = *(tightPhotonFlag.product());
236 
237  // Get the reconstructed muons
238  bool validMuons=true;
239  Handle<reco::MuonCollection> muonHandle;
240  reco::MuonCollection muonCollection;
241  e.getByLabel(muonProducer_, muonCollection_ , muonHandle);
242  if ( !muonHandle.isValid()) {
243  edm::LogInfo("ZToMuMuGammaAnalyzer") << "Error! Can't get the product "<< muonCollection_ << endl;
244  validMuons=false;
245  }
246  if(validMuons) muonCollection = *(muonHandle.product());
247 
248  // Get the beam spot
250  e.getByLabel("offlineBeamSpot", bsHandle);
251  if (!bsHandle.isValid()) {
252  edm::LogError("TrackerOnlyConversionProducer") << "Error! Can't get the product primary Vertex Collection "<< "\n";
253  return;
254  }
255  const reco::BeamSpot &thebs = *bsHandle.product();
256 
257 
258 
259 
260  //Prepare list of photon-related HLT filter names
261  vector<int> Keys;
262  for(uint filterIndex=0;filterIndex<triggerEvent.sizeFilters();++filterIndex){ //loop over all trigger filters in event (i.e. filters passed)
263  string label = triggerEvent.filterTag(filterIndex).label();
264  if(label.find( "Photon" ) != string::npos ) { //get photon-related filters
265  for(uint filterKeyIndex=0;filterKeyIndex<triggerEvent.filterKeys(filterIndex).size();++filterKeyIndex){ //loop over keys to objects passing this filter
266  Keys.push_back(triggerEvent.filterKeys(filterIndex)[filterKeyIndex]); //add keys to a vector for later reference
267  }
268  }
269  }
270 
271  // sort Keys vector in ascending order
272  // and erases duplicate entries from the vector
273  sort(Keys.begin(),Keys.end());
274  for ( uint i=0 ; i<Keys.size() ; )
275  {
276  if (i!=(Keys.size()-1))
277  {
278  if (Keys[i]==Keys[i+1]) Keys.erase(Keys.begin()+i+1) ;
279  else ++i ;
280  }
281  else ++i ;
282  }
283 
284 
285 
287  if ( muonCollection.size() < 2 ) return;
288 
289  for( reco::MuonCollection::const_iterator iMu = muonCollection.begin(); iMu != muonCollection.end(); iMu++) {
290  if ( !basicMuonSelection (*iMu) ) continue;
291 
292  for( reco::MuonCollection::const_iterator iMu2 = iMu+1; iMu2 != muonCollection.end(); iMu2++) {
293  if ( !basicMuonSelection (*iMu2) ) continue;
294  if ( iMu->charge()*iMu2->charge() > 0) continue;
295 
296  if ( !muonSelection(*iMu,thebs) && !muonSelection(*iMu2,thebs) ) continue;
297 
298 
299 
300  float mumuMass = mumuInvMass(*iMu,*iMu2) ;
301  if ( mumuMass < minMumuInvMass_ || mumuMass > maxMumuInvMass_ ) continue;
302 
303  h1_mumuInvMass_ -> Fill (mumuMass);
304 
305  if ( photonCollection.size() < 1 ) continue;
306 
307  reco::Muon nearMuon;
308  reco::Muon farMuon;
309  for( reco::PhotonCollection::const_iterator iPho = photonCollection.begin(); iPho != photonCollection.end(); iPho++) {
310  if ( !photonSelection (*iPho) ) continue;
311 
312 
314  double dr1 = deltaR(*iMu, *iPho);
315  double dr2 = deltaR(*iMu2,*iPho);
316  double drNear = dr1;
317  if (dr1 < dr2) {
318  nearMuon =*iMu ; farMuon = *iMu2; drNear = dr1;
319  } else {
320  nearMuon = *iMu2; farMuon = *iMu; drNear = dr2;
321  }
322 
323  if ( nearMuon.isolationR03().hadEt > nearMuonHcalIso_ ) continue;
324  if ( farMuon.isolationR03().sumPt > farMuonTrackIso_ ) continue;
325  if ( farMuon.isolationR03().emEt > farMuonEcalIso_ ) continue;
326  if ( farMuon.pt() < farMuonMinPt_ ) continue;
327  if ( drNear > nearMuonDr_) continue;
328 
329 
330  float mumuGammaMass = mumuGammaInvMass(*iMu,*iMu2,*iPho) ;
331  if ( mumuGammaMass < minMumuGammaInvMass_ || mumuGammaMass > maxMumuGammaInvMass_ ) continue;
332 
333  h1_mumuGammaInvMass_ ->Fill (mumuGammaMass);
334 
335 
336  }
337 
338  }
339 
340  }
341 
342 
343 }//End of Analyze method
344 
346 {
347  if(!standAlone_){
348 
349  dbe_->setCurrentFolder("Egamma/PhotonAnalyzer/ZToMuMuGamma");
350 
351 
352  }
353 
354 }
355 
356 
358 {
359  //dbe_->showDirStructure();
360  if(standAlone_){
361  dbe_->setCurrentFolder("Egamma/PhotonAnalyzer/ZToMuMuGamma");
362 
363  dbe_->save(outputFileName_);
364  }
365 
366 
367 }
368 
370  bool result=true;
371  if (!mu.innerTrack().isNonnull()) result=false;
372  if (!mu.globalTrack().isNonnull()) result=false;
373  if ( !mu.isGlobalMuon() ) result=false;
374  if ( mu.pt() < muonMinPt_ ) result=false;
375  if ( fabs(mu.eta())>2.4 ) result=false;
376 
377  int pixHits=0;
378  int tkHits=0;
379  if ( mu.innerTrack().isNonnull() ) {
380  pixHits=mu.innerTrack()->hitPattern().numberOfValidPixelHits();
381  tkHits=mu.innerTrack()->hitPattern().numberOfValidStripHits();
382  }
383 
384  if ( pixHits+tkHits < minPixStripHits_ ) result=false;
385 
386 
387  return result;
388 }
389 
391  bool result=true;
392  if ( mu.globalTrack()->normalizedChi2() > muonMaxChi2_ ) result=false;
393  if ( fabs( mu.globalTrack()->dxy(beamSpot)) > muonMaxDxy_ ) result=false;
394  if ( mu.numberOfMatches() < muonMatches_ ) result=false;
395 
396  if ( mu.track()-> hitPattern().numberOfValidPixelHits() < validPixHits_ ) result=false;
397  if ( mu.globalTrack()->hitPattern().numberOfValidMuonHits() < validMuonHits_ ) result=false;
398  if ( !mu.isTrackerMuon() ) result=false;
399  // track isolation
400  if ( mu.isolationR03().sumPt > muonTrackIso_ ) result=false;
401  if ( fabs(mu.eta())> muonTightEta_ ) result=false;
402 
403 
404  return result;
405 }
406 
407 
409  bool result=true;
410  if ( pho.pt() < photonMinEt_ ) result=false;
411  if ( fabs(pho.eta())> photonMaxEta_ ) result=false;
412  if ( pho.isEBEEGap() ) result=false;
413  // if ( pho.trkSumPtHollowConeDR04() > photonTrackIso_ ) result=false; // check how to exclude the muon track (which muon track).
414 
415 
416  return result;
417 }
418 
419 
420 
422  {
423  math::XYZTLorentzVector p12 = mu1.p4()+mu2.p4() ;
424  float mumuMass2 = p12.Dot(p12) ;
425  float invMass = sqrt(mumuMass2) ;
426  return invMass ;
427  }
428 
430  {
431  math::XYZTLorentzVector p12 = mu1.p4()+mu2.p4()+pho.p4() ;
432  float Mass2 = p12.Dot(p12) ;
433  float invMass = sqrt(Mass2) ;
434  return invMass ;
435  }
436 
float hadEt
hcal sum-Et
Definition: MuonIsolation.h:9
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
float sumPt
sum-pt of tracks
Definition: MuonIsolation.h:7
bool muonSelection(const reco::Muon &m, const reco::BeamSpot &bs)
The single EDProduct to be saved for each event (AOD case)
Definition: TriggerEvent.h:27
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:717
Definition: deltaR.h:51
virtual TrackRef innerTrack() const
Definition: Muon.h:49
bool isTrackerMuon() const
Definition: Muon.h:212
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE")
Definition: DQMStore.cc:2113
bool isEBEEGap() const
true if photon is in boundary between EB and EE
Definition: Photon.h:130
virtual TrackRef track() const
reference to a Track
Definition: Muon.h:50
bool isGlobalMuon() const
Definition: Muon.h:211
virtual double eta() const
momentum pseudorapidity
std::vector< Muon > MuonCollection
collection of Muon objects
Definition: MuonFwd.h:9
virtual void analyze(const edm::Event &, const edm::EventSetup &)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:30
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
const LorentzVector & p4(P4type type) const
Definition: Photon.cc:197
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
T sqrt(T t)
Definition: SSEVec.h:46
float emEt
ecal sum-Et
Definition: MuonIsolation.h:8
tuple result
Definition: query.py:137
float mumuInvMass(const reco::Muon &m1, const reco::Muon &m2)
const int mu
Definition: Constants.h:23
ZToMuMuGammaAnalyzer(const edm::ParameterSet &)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
DQMStore * dbe_
int numberOfMatches(ArbitrationType type=SegmentAndTrackArbitration) const
get number of chambers with matched segments
Definition: Muon.cc:56
bool etMin(const PFCandidate &cand, double cut)
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
bool photonSelection(const reco::Photon &p)
virtual double pt() const
transverse momentum
std::vector< size_type > Keys
std::vector< Photon > PhotonCollection
collectin of Photon objects
Definition: PhotonFwd.h:9
edm::EventID id() const
Definition: EventBase.h:56
bool basicMuonSelection(const reco::Muon &m)
virtual void endRun(const edm::Run &, const edm::EventSetup &)
virtual const LorentzVector & p4() const
four-momentum Lorentz vector
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
float mumuGammaInvMass(const reco::Muon &mu1, const reco::Muon &mu2, const reco::Photon &pho)
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:429
const MuonIsolation & isolationR03() const
Definition: Muon.h:159
Definition: Run.h:33
virtual TrackRef globalTrack() const
reference to Track reconstructed in both tracked and muon detector
Definition: Muon.h:55