CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ZMuMuIsolationAnalyzer.cc
Go to the documentation of this file.
22 #include "TH1.h"
23 #include "TH2.h"
24 #include "TMath.h"
25 #include <vector>
26 #include <string>
27 #include <iostream>
28 #include <sstream>
29 
30 using namespace edm;
31 using namespace std;
32 using namespace reco;
33 using namespace isodeposit;
34 
36 public:
38 private:
39  virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
40  virtual void endJob() override;
42  double dRVeto;
43  double dRTrk, dREcal, dRHcal;
44  double ptThreshold, etEcalThreshold, etHcalThreshold;
45  double alpha, beta;
46  double pt, eta;
47  double iso_cut;
48 
49  TH1F * h_IsoZ_tk,* h_IsoW_tk,* h_IsoOther_tk ;
50  TH1F * h_IsoZ_ecal,* h_IsoW_ecal,* h_IsoOther_ecal ;
51  TH1F * h_IsoZ_hcal,* h_IsoW_hcal,* h_IsoOther_hcal ;
52  TH1F * IsoZ,* IsoW,* IsoOther ;
53  TH1F * TkrPt,* EcalEt,* HcalEt ;
54  TH1F * EcalEtZ, * HcalEtZ;
55 
56  TH1F * Z_eta,* W_eta,* Other_eta;
57  TH1F * Z_eta_postSelection,* W_eta_postSelection,* Other_eta_postSelection;
58  TH1F * Z_pt,* W_pt,* Other_pt;
59  TH1F * Z_pt_postSelection,* W_pt_postSelection,* Other_pt_postSelection;
60 
61  enum MuTag { muFromZ, muFromW, muFromOther };
62  template<typename T>
63  MuTag muTag(const T & mu) const;
64  void Deposits(const pat::IsoDeposit* isodep, double dR_max, TH1F* hist);
65  void histo(TH1F* hist, const char* cx, const char* cy) const;
66 };
67 
68 template<typename T>
70  GenParticleRef p = mu.genParticleRef();
71  if(p.isNull()){
72  // cout<<"genParticleRef is null "<<endl;
73  return muFromOther;
74 }
75  int sizem = p->numberOfMothers();
76  if(sizem != 1) {
77  //cout<<"number of mothers !=1 "<<endl;
78  return muFromOther;
79  }
80  const Candidate * moth1 = p->mother();
81  if(moth1 == 0) {
82  return muFromOther;
83  //cout<<"no mother "<<endl;
84  }
85  int pdgId1 = moth1->pdgId();
86  if(abs(pdgId1)!=13){
87  return muFromOther;
88  //cout<<"mother is not a muon"<<endl;
89  }
90  const Candidate * moth2 = moth1->mother();
91  if(moth2 == 0) {
92  return muFromOther;
93  //cout<<"no mother "<<endl;
94 }
95  int pdgId2 = moth2->pdgId();
96  if(pdgId2 == 23) {
97  //cout<<" muon from Z"<<endl;
98  return muFromZ;
99  }
100  if(abs(pdgId2)==24) return muFromW;
101  else {
102  //cout<<" muon from other"<<endl;
103  return muFromOther;
104  }
105 }
106 
107 void ZMuMuIsolationAnalyzer::Deposits(const pat::IsoDeposit* isodep,double dR_max,TH1F* hist){
108  for(IsoDeposit::const_iterator it= isodep->begin(); it!= isodep->end(); ++it){
109  if(it->dR()<dR_max) {
110  double theta= 2*(TMath::ATan(TMath::Exp(-(it->eta() ) ) ) );
111  // double theta= 2;
112  hist->Fill(it->value()/TMath::Sin(theta));
113  }
114  }
115 }
116 
117 void ZMuMuIsolationAnalyzer::histo(TH1F* hist, const char* cx, const char*cy) const{
118  hist->GetXaxis()->SetTitle(cx);
119  hist->GetYaxis()->SetTitle(cy);
120  hist->GetXaxis()->SetTitleOffset(1);
121  hist->GetYaxis()->SetTitleOffset(1.2);
122  hist->GetXaxis()->SetTitleSize(0.04);
123  hist->GetYaxis()->SetTitleSize(0.04);
124  hist->GetXaxis()->SetLabelSize(0.03);
125  hist->GetYaxis()->SetLabelSize(0.03);
126 }
127 
129  srcToken(consumes<CandidateView>(pset.getParameter<InputTag>("src"))),
130  dRVeto(pset.getUntrackedParameter<double>("veto")),
131  dRTrk(pset.getUntrackedParameter<double>("deltaRTrk")),
132  dREcal(pset.getUntrackedParameter<double>("deltaREcal")),
133  dRHcal(pset.getUntrackedParameter<double>("deltaRHcal")),
134  ptThreshold(pset.getUntrackedParameter<double>("ptThreshold")),
135  etEcalThreshold(pset.getUntrackedParameter<double>("etEcalThreshold")),
136  etHcalThreshold(pset.getUntrackedParameter<double>("etHcalThreshold")),
137  alpha(pset.getUntrackedParameter<double>("alpha")),
138  beta(pset.getUntrackedParameter<double>("beta")),
139  pt(pset.getUntrackedParameter<double>("pt")),
140  eta(pset.getUntrackedParameter<double>("eta")),
141  iso_cut(pset.getUntrackedParameter<double>("isoCut")) {
143  std::ostringstream str1,str2,str3,str4,str5,str6,str7,str8,str9,str10,n_tracks;
144  str1 << "muons from Z with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk;
145  str2 << "muons from W with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk;
146  str3 << "muons from Others with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk;
147  str4 << "muons from Z with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<" (alpha = "<<alpha<<" , "<<"beta = "<<beta<<" )";
148  str5 << "muons from W with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<" (alpha = "<<alpha<<" , "<<"beta = "<<beta<<" )";
149  str6 << "muons from Other with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<" (alpha = "<<alpha<<" , "<<"beta = "<<beta<<" )";
150  n_tracks <<"Number of tracks for muon with p_{t} > " << ptThreshold <<" and #Delta R < "<<dRTrk<< " GeV/c";
151  str7<<"Isolation Vs p_{t} with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<"(Tracker)";
152  str8<<"Isolation Vs p_{t} with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<"(Ecal)";
153  str9<<"Isolation Vs p_{t} with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<"(Hcal)";
154  str10<<"Isolation Vs p_{t} with p_{t} > " << ptThreshold << " GeV/c"<<" and #Delta R < "<<dRTrk<<" (alpha = "<<alpha<<" , "<<"beta = "<<beta<<" )";
155  h_IsoZ_tk = fs->make<TH1F>("ZIso_Tk",str1.str().c_str(),100,0.,20.);
156  h_IsoW_tk = fs->make<TH1F>("WIso_Tk",str2.str().c_str(),100,0.,20.);
157  h_IsoOther_tk = fs->make<TH1F>("otherIso_Tk",str3.str().c_str(),100,0.,20.);
158  h_IsoZ_ecal = fs->make<TH1F>("ZIso_ecal",str1.str().c_str(),100,0.,20.);
159  h_IsoW_ecal = fs->make<TH1F>("WIso_ecal",str2.str().c_str(),100,0.,20.);
160  h_IsoOther_ecal = fs->make<TH1F>("otherIso_ecal",str3.str().c_str(),100,0.,20.);
161  h_IsoZ_hcal = fs->make<TH1F>("ZIso_hcal",str1.str().c_str(),100,0.,20.);
162  h_IsoW_hcal = fs->make<TH1F>("WIso_hcal",str2.str().c_str(),100,0.,20.);
163  h_IsoOther_hcal = fs->make<TH1F>("otherIso_hcal",str3.str().c_str(),100,0.,20.);
164  IsoZ = fs->make<TH1F>("ZIso",str4.str().c_str(),100,0.,20.);
165  IsoW = fs->make<TH1F>("WIso",str5.str().c_str(),100,0.,20.);
166  IsoOther = fs->make<TH1F>("otherIso",str6.str().c_str(),100,0.,20.);
167 
168 
169  Z_eta = fs->make<TH1F>("Z_eta","#eta distribution for muons coming from Z",40,-eta,eta);
170  W_eta = fs->make<TH1F>("W_eta","#eta distribution for muons coming from W",40,-eta,eta);
171  Other_eta = fs->make<TH1F>("Other_eta","#eta distribution for muons coming from other",40,-eta,eta);
172  Z_eta_postSelection = fs->make<TH1F>("Z_eta_postSelection","#eta distribution for muons coming from Z after iso selection",40,-eta,eta);
173  W_eta_postSelection = fs->make<TH1F>("W_eta_postSelection","#eta distribution for muons coming from W after iso selection",40,-eta,eta);
174  Other_eta_postSelection = fs->make<TH1F>("Other_eta_postSelection","#eta distribution for muons coming from other after iso selection",40,-eta,eta);
175 
176  Z_pt = fs->make<TH1F>("Z_pt","p_{T} distribution for muons coming from Z",40,pt,150.);
177  W_pt = fs->make<TH1F>("W_pt","p_{T} distribution for muons coming from W",40,pt,150.);
178  Other_pt = fs->make<TH1F>("Other_pt","p_{T} distribution for muons coming from other",40,pt,150.);
179  Z_pt_postSelection = fs->make<TH1F>("Z_pt_postSelection","p_{T} distribution for muons coming from Z after iso selection",40,pt,150.);
180  W_pt_postSelection = fs->make<TH1F>("W_pt_postSelection","p_{t} distribution for muons coming from W after iso selection",40,pt,150.);
181  Other_pt_postSelection = fs->make<TH1F>("Other_pt_postSelection","p_{t} distribution for muons coming from other after iso selection",40,pt,150.);
182 
183 
184  TkrPt = fs->make<TH1F>("TkrPt","IsoDeposit p distribution in the Tracker",100,0.,10.);
185  EcalEt = fs->make<TH1F>("EcalEt","IsoDeposit E distribution in the Ecal",100,0.,5.);
186  HcalEt = fs->make<TH1F>("HcalEt","IsoDeposit E distribution in the Hcal",100,0.,5.);
187 
188  EcalEtZ = fs->make<TH1F>("VetoEcalEt"," #Sigma E_{T} deposited in veto cone in the Ecal",100,0.,10.);
189  HcalEtZ = fs->make<TH1F>("VetoHcalEt"," #Sigma E_{T} deposited in veto cone in the Hcal",100,0.,10.);
190 }
191 
194  event.getByToken(srcToken,dimuons);
195 
196  for(unsigned int i=0; i < dimuons->size(); ++ i ) {
197  const Candidate & zmm = (* dimuons)[i];
198  const Candidate * dau0 = zmm.daughter(0);
199  const Candidate * dau1 = zmm.daughter(1);
200  const pat::Muon & mu0 = dynamic_cast<const pat::Muon &>(*dau0->masterClone());
201  const pat::GenericParticle & mu1 = dynamic_cast<const pat::GenericParticle &>(*dau1->masterClone());
202 
203  const pat::IsoDeposit * muTrackIso =mu0.isoDeposit(pat::TrackIso);
204  const pat::IsoDeposit * tkTrackIso =mu1.isoDeposit(pat::TrackIso);
205  const pat::IsoDeposit * muEcalIso =mu0.isoDeposit(pat::EcalIso);
206  const pat::IsoDeposit * tkEcalIso =mu1.isoDeposit(pat::EcalIso);
207  const pat::IsoDeposit * muHcalIso =mu0.isoDeposit(pat::HcalIso);
208  const pat::IsoDeposit * tkHcalIso =mu1.isoDeposit(pat::HcalIso);
209 
210 
211  if(mu0.pt() > pt && mu1.pt() > pt && abs(mu0.eta()) < eta && abs(mu1.eta()) < eta){
212 
213  Direction muDir = Direction(mu0.eta(),mu0.phi());
214  Direction tkDir = Direction(mu1.eta(),mu1.phi());
215 
216  IsoDeposit::AbsVetos vetos_mu;
217  vetos_mu.push_back(new ConeVeto( muDir, dRVeto ));
218  vetos_mu.push_back(new ThresholdVeto( ptThreshold ));
219 
221  vetos_tk.push_back(new ConeVeto( tkDir, dRVeto ));
222  vetos_tk.push_back(new ThresholdVeto( ptThreshold ));
223 
224  reco::IsoDeposit::AbsVetos vetos_mu_ecal;
225  vetos_mu_ecal.push_back(new ConeVeto( muDir, 0. ));
226  vetos_mu_ecal.push_back(new ThresholdVeto( etEcalThreshold ));
227 
228  reco::IsoDeposit::AbsVetos vetos_tk_ecal;
229  vetos_tk_ecal.push_back(new ConeVeto( tkDir, 0. ));
230  vetos_tk_ecal.push_back(new ThresholdVeto( etEcalThreshold ));
231 
232  reco::IsoDeposit::AbsVetos vetos_mu_hcal;
233  vetos_mu_hcal.push_back(new ConeVeto( muDir, 0. ));
234  vetos_mu_hcal.push_back(new ThresholdVeto( etHcalThreshold ));
235 
236  reco::IsoDeposit::AbsVetos vetos_tk_hcal;
237  vetos_tk_hcal.push_back(new ConeVeto( tkDir, 0. ));
238  vetos_tk_hcal.push_back(new ThresholdVeto( etHcalThreshold ));
239  MuTag tag_mu = muTag(mu0);
240  MuTag tag_track = muTag(mu1);
241 
242  double Tk_isovalue = TMath::Max(muTrackIso->sumWithin(dRTrk,vetos_mu),tkTrackIso->sumWithin(dRTrk, vetos_tk));
243  double Ecal_isovalue = TMath::Max(muEcalIso->sumWithin(dREcal,vetos_mu_ecal),tkEcalIso->sumWithin(dREcal, vetos_tk_ecal));
244  double Hcal_isovalue = TMath::Max(muHcalIso->sumWithin(dRHcal,vetos_mu_hcal),tkHcalIso->sumWithin(dRHcal, vetos_tk_hcal));
245  EcalEtZ->Fill(muEcalIso->candEnergy());
246  EcalEtZ->Fill(tkEcalIso->candEnergy());
247  HcalEtZ->Fill(muHcalIso->candEnergy());
248  HcalEtZ->Fill(tkHcalIso->candEnergy());
249 
250  double iso_value0 = alpha*((0.5*(1+beta)* muEcalIso->sumWithin(dREcal,vetos_mu_ecal) ) + (0.5*(1-beta)*muHcalIso->sumWithin(dRHcal,vetos_mu_hcal) ) ) +(1-alpha)*muTrackIso->sumWithin(dRTrk,vetos_mu);
251  double iso_value1 = alpha*((0.5*(1+beta)* tkEcalIso->sumWithin(dREcal,vetos_tk_ecal) ) + (0.5*(1-beta)*tkHcalIso->sumWithin(dRHcal,vetos_tk_hcal) ) ) +(1-alpha)*tkTrackIso->sumWithin(dRTrk,vetos_tk);
252 
253  double iso_value=TMath::Max(iso_value0,iso_value1);
254 
255  if(tag_mu==muFromZ && tag_track==muFromZ){
256  h_IsoZ_tk->Fill(Tk_isovalue);
257  h_IsoZ_ecal->Fill(Ecal_isovalue);
258  h_IsoZ_hcal->Fill(Hcal_isovalue);
259  IsoZ->Fill(iso_value);
260 
261  Z_eta->Fill(mu0.eta());
262  Z_eta->Fill(mu1.eta());
263  Z_pt->Fill(mu0.pt());
264  Z_pt->Fill(mu1.pt());
265 
266  if(iso_value0<iso_cut) {
267  Z_pt_postSelection->Fill(mu0.pt());
268  Z_eta_postSelection->Fill(mu0.eta());
269  }
270  if(iso_value1<iso_cut){
271  Z_pt_postSelection->Fill(mu1.pt());
272  Z_eta_postSelection->Fill(mu1.eta());
273  }
274 
275  Deposits(muTrackIso,dRTrk,TkrPt);
276  Deposits(muEcalIso,dREcal,EcalEt);
277  Deposits(muHcalIso,dRHcal,HcalEt);
278  Deposits(tkTrackIso,dRTrk,TkrPt);
279  Deposits(tkEcalIso,dREcal,EcalEt);
280  Deposits(tkHcalIso,dRHcal,HcalEt);
281  }
282  if(tag_mu==muFromW || tag_track==muFromW){
283  h_IsoW_tk->Fill(Tk_isovalue);
284  h_IsoW_ecal->Fill(Ecal_isovalue);
285  h_IsoW_hcal->Fill(Hcal_isovalue);
286  IsoW->Fill(iso_value);
287 
288  W_eta->Fill(mu0.eta());
289  W_eta->Fill(mu1.eta());
290  W_pt->Fill(mu0.pt());
291  W_pt->Fill(mu1.pt());
292 
293  if(iso_value0<iso_cut) {
294  W_pt_postSelection->Fill(mu0.pt());
295  W_eta_postSelection->Fill(mu0.eta());
296  }
297  if(iso_value1<iso_cut) {
298  W_pt_postSelection->Fill(mu1.pt());
299  W_eta_postSelection->Fill(mu1.eta());
300  }
301 
302  Deposits(muTrackIso,dRTrk,TkrPt);
303  Deposits(muEcalIso,dREcal,EcalEt);
304  Deposits(muHcalIso,dRHcal,HcalEt);
305  Deposits(tkTrackIso,dRTrk,TkrPt);
306  Deposits(tkEcalIso,dREcal,EcalEt);
307  Deposits(tkHcalIso,dRHcal,HcalEt);
308  }
309  else{
310  h_IsoOther_tk->Fill(Tk_isovalue);
311  h_IsoOther_ecal->Fill(Ecal_isovalue);
312  h_IsoOther_hcal->Fill(Hcal_isovalue);
313  IsoOther->Fill(iso_value);
314 
315  Other_eta->Fill(mu0.eta());
316  Other_eta->Fill(mu1.eta());
317  Other_pt->Fill(mu0.pt());
318  Other_pt->Fill(mu1.pt());
319 
320  if(iso_value0<iso_cut) {
321  Other_pt_postSelection->Fill(mu0.pt());
322  Other_eta_postSelection->Fill(mu0.eta());
323  }
324  if(iso_value1<iso_cut) {
325  Other_pt_postSelection->Fill(mu1.pt());
326  Other_eta_postSelection->Fill(mu1.eta());
327  }
328 
329  Deposits(muTrackIso,dRTrk,TkrPt);
330  Deposits(muEcalIso,dREcal,EcalEt);
331  Deposits(muHcalIso,dRHcal,HcalEt);
332  Deposits(tkTrackIso,dRTrk,TkrPt);
333  Deposits(tkEcalIso,dREcal,EcalEt);
334  Deposits(tkHcalIso,dRHcal,HcalEt);
335  }
336  }
337  }
338 
339  histo(h_IsoZ_tk,"#Sigma p_{T}","Events");
340  histo(h_IsoW_tk,"#Sigma p_{T}","Events");
341  histo(h_IsoOther_tk,"#Sigma p_{T}","#Events");
342  histo(h_IsoZ_ecal,"#Sigma E_{t}","Events");
343  histo(h_IsoW_ecal,"#Sigma E_{t}","Events");
344  histo(h_IsoOther_ecal,"#Sigma E_{t}","Events");
345  histo(h_IsoZ_hcal,"#Sigma E_{t}","Events");
346  histo(h_IsoW_hcal,"#Sigma E_{t}","Events");
347  histo(h_IsoOther_hcal,"#Sigma E_{t}","Events");
348  histo(TkrPt,"p ","");
349  histo(EcalEt,"E ","");
350  histo(HcalEt,"E ","");
351  histo(HcalEtZ,"E_{T}","");
352  histo(EcalEtZ,"E_{T}","");
353 }
354 
356 }
357 
359 
const double beta
double candEnergy() const
Get energy or pT attached to cand trajectory.
Definition: IsoDeposit.h:136
int i
Definition: DBlmapReader.cc:9
float alpha
Definition: AMPTWrapper.h:95
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
EDGetTokenT< CandidateView > srcToken
Geom::Theta< T > theta() const
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
const IsoDeposit * isoDeposit(IsolationKeys key) const
Returns the IsoDeposit associated with some key, or a null pointer if it is not available.
virtual void endJob() override
T eta() const
const IsoDeposit * isoDeposit(IsolationKeys key) const
Returns the IsoDeposit associated with some key, or a null pointer if it is not available.
Definition: Lepton.h:165
virtual double eta() const
momentum pseudorapidity
virtual double pt() const
transverse momentum
double sumWithin(double coneSize, const AbsVetos &vetos=AbsVetos(), bool skipDepositVeto=false) const
Definition: IsoDeposit.cc:138
Analysis-level Generic Particle class (e.g. for hadron or muon not fully reconstructed) ...
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ZMuMuIsolationAnalyzer(const edm::ParameterSet &pset)
const int mu
Definition: Constants.h:22
void Deposits(const pat::IsoDeposit *isodep, double dR_max, TH1F *hist)
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
bool isNull() const
Checks for null.
Definition: Ref.h:247
const_iterator begin() const
Definition: IsoDeposit.h:163
virtual int pdgId() const =0
PDG identifier.
T Max(T a, T b)
Definition: MathUtil.h:44
std::vector< AbsVeto * > AbsVetos
Definition: IsoDeposit.h:40
MuTag muTag(const T &mu) const
isodeposit::AbsVetos AbsVetos
Definition: IsoDeposit.h:51
const_iterator end() const
Definition: IsoDeposit.h:164
long double T
void histo(TH1F *hist, const char *cx, const char *cy) const
virtual double phi() const
momentum azimuthal angle
Analysis-level muon class.
Definition: Muon.h:49
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
virtual void analyze(const edm::Event &event, const edm::EventSetup &setup) override
virtual const CandidateBaseRef & masterClone() const =0