CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ZHistogrammer.cc
Go to the documentation of this file.
1 #include "FWCore/Framework/interface/EDAnalyzer.h"
2 #include "FWCore/Utilities/interface/InputTag.h"
3 #include "TH1.h"
4 
5 class ZHistogrammer : public edm::EDAnalyzer {
6 public:
7  ZHistogrammer(const edm::ParameterSet& pset);
8 private:
9  virtual void analyze(const edm::Event& event, const edm::EventSetup& setup);
10  edm::InputTag z_, gen_, match_;
11  unsigned int nbinsMass_, nbinsPt_, nbinsAng_, nbinsMassRes_;
12  double massMax_, ptMax_, angMax_, massResMax_;
13  TH1F *h_nZ_, *h_mZ_, *h_ptZ_, *h_phiZ_, *h_thetaZ_, *h_etaZ_, *h_rapidityZ_;
14  TH1F *h_invmMuMu_;
15  TH1F *h_nZMC_, *h_mZMC_, *h_ptZMC_, *h_phiZMC_, *h_thetaZMC_, *h_etaZMC_, *h_rapidityZMC_;
16  TH1F *h_invmMuMuMC_;
17  //TH1F *h_mZ2vs3MC_, *h_ptZ2vs3MC_, *h_phiZ2vs3MC_, *h_thetaZ2vs3MC_, *h_etaZ2vs3MC_, *h_rapidityZ2vs3MC_;
18  TH1F *h_mResZ_, *h_ptResZ_, *h_phiResZ_, *h_thetaResZ_, *h_etaResZ_, *h_rapidityResZ_;
19  TH1F *h_mResZMuMu_, *h_mRatioZMuMu_;
20  TH1F *h_mResZMuMuMC_, *h_mRatioZMuMuMC_;
21 };
22 
23 #include "DataFormats/Candidate/interface/Candidate.h"
24 #include "DataFormats/Candidate/interface/CandidateFwd.h"
25 #include "DataFormats/Candidate/interface/CandMatchMap.h"
26 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
27 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
28 #include "FWCore/Framework/interface/Event.h"
29 #include "DataFormats/Common/interface/Handle.h"
30 #include "FWCore/ParameterSet/interface/ParameterSet.h"
31 #include "FWCore/ServiceRegistry/interface/Service.h"
32 #include "CommonTools/UtilAlgos/interface/TFileService.h"
33 #include <cmath>
34 #include <iostream>
35 
36 using namespace std;
37 using namespace reco;
38 using namespace edm;
39 
40 ZHistogrammer::ZHistogrammer(const ParameterSet& pset) :
41  z_(pset.getParameter<InputTag>("z")),
42  gen_(pset.getParameter<InputTag>("gen")),
43  match_(pset.getParameter<InputTag>("match")),
44  nbinsMass_(pset.getUntrackedParameter<unsigned int>("nbinsMass")),
45  nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")),
46  nbinsAng_(pset.getUntrackedParameter<unsigned int>("nbinsAng")),
47  nbinsMassRes_(pset.getUntrackedParameter<unsigned int>("nbinsMassRes")),
48  massMax_(pset.getUntrackedParameter<double>("massMax")),
49  ptMax_(pset.getUntrackedParameter<double>("ptMax")),
50  angMax_(pset.getUntrackedParameter<double>("angMax")),
51  massResMax_(pset.getUntrackedParameter<double>("massResMax")) {
52  cout << ">>> Z Histogrammer constructor" << endl;
53  Service<TFileService> fs;
54  TFileDirectory ZHisto = fs->mkdir( "ZRecoHisto" );
55  TFileDirectory ZMCHisto = fs->mkdir( "ZMCHisto" );
56  TFileDirectory ZResHisto = fs->mkdir( "ZResolutionHisto" );
57  //TFileDirectory Z2vs3MCHisto = fs->mkdir( "Z2vs3MCHisto" );
58  h_nZ_ = ZHisto.make<TH1F>("ZNumber", "number of Z particles", 11, -0.5, 10.5);
59  h_mZ_ = ZHisto.make<TH1F>("ZMass", "Z mass (GeV/c^{2})", nbinsMass_, 0, massMax_);
60  h_ptZ_ = ZHisto.make<TH1F>("ZPt", "Z p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
61  h_phiZ_ = ZHisto.make<TH1F>("ZPhi", "Z #phi", nbinsAng_, -angMax_, angMax_);
62  h_thetaZ_ = ZHisto.make<TH1F>("Ztheta", "Z #theta", nbinsAng_, 0, angMax_);
63  h_etaZ_ = ZHisto.make<TH1F>("ZEta", "Z #eta", nbinsAng_, -angMax_, angMax_);
64  h_rapidityZ_ = ZHisto.make<TH1F>("ZRapidity", "Z rapidity", nbinsAng_, -angMax_, angMax_);
65  h_invmMuMu_ = ZHisto.make<TH1F>("MuMuMass", "#mu #mu invariant mass",
66  nbinsMass_, 0, massMax_);
67  h_nZMC_ = ZMCHisto.make<TH1F>("ZMCNumber", "number of Z MC particles", 11, -0.5, 10.5);
68  h_mZMC_ = ZMCHisto.make<TH1F>("ZMCMass", "Z MC mass (GeV/c^{2})", nbinsMass_, 0, massMax_);
69  h_ptZMC_ = ZMCHisto.make<TH1F>("ZMCPt", "Z MC p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
70  h_phiZMC_ = ZMCHisto.make<TH1F>("ZMCPhi", "Z MC #phi", nbinsAng_, -angMax_, angMax_);
71  h_thetaZMC_ = ZMCHisto.make<TH1F>("ZMCTheta", "Z MC #theta", nbinsAng_, 0, angMax_);
72  h_etaZMC_ = ZMCHisto.make<TH1F>("ZMCEta", "Z MC #eta", nbinsAng_, -angMax_, angMax_);
73  h_rapidityZMC_ = ZMCHisto.make<TH1F>("ZMCRapidity", "Z MC rapidity",
74  nbinsAng_, -angMax_, angMax_);
75  h_invmMuMuMC_ = ZMCHisto.make<TH1F>("MuMuMCMass", "#mu #mu MC invariant mass",
76  nbinsMass_, 0, massMax_);
77  /*
78  h_mZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCMass", "Z MC st 2 vs st 3 mass (GeV/c^{2})",
79  nbinsMassRes_, -massResMax_, massResMax_);
80  h_ptZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCPt", "Z MC st 2 vs st 3 p_{t} (GeV/c)",
81  nbinsPt_, -ptMax_, ptMax_);
82  h_phiZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCPhi", "Z MC st 2 vs st 3 #phi",
83  nbinsAng_, -angMax_, angMax_);
84  h_thetaZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCTheta", "Z MC st 2 vs st 3 #theta",
85  nbinsAng_, -angMax_, angMax_);
86  h_etaZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCEta", "Z MC st 2 vs st 3 #eta",
87  nbinsAng_, -angMax_, angMax_);
88  h_rapidityZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCRapidity", "Z MC st 2 vs st 3 rapidity",
89  nbinsAng_, -angMax_, angMax_);
90  */
91  h_mResZ_ = ZResHisto.make<TH1F>("ZMassResolution", "Z mass Resolution (GeV/c^{2})",
92  nbinsMassRes_, -massResMax_, massResMax_);
93  h_ptResZ_ = ZResHisto.make<TH1F>("ZPtResolution", "Z p_{t} Resolution (GeV/c)",
94  nbinsPt_, -ptMax_, ptMax_);
95  h_phiResZ_ = ZResHisto.make<TH1F>("ZPhiResolution", "Z #phi Resolution",
96  nbinsAng_, -angMax_, angMax_);
97  h_thetaResZ_ = ZResHisto.make<TH1F>("ZThetaResolution", "Z #theta Resolution",
98  nbinsAng_, -angMax_, angMax_);
99  h_etaResZ_ = ZResHisto.make<TH1F>("ZEtaResolution", "Z #eta Resolution",
100  nbinsAng_, -angMax_, angMax_);
101  h_rapidityResZ_ = ZResHisto.make<TH1F>("ZRapidityResolution", "Z rapidity Resolution",
102  nbinsAng_, -angMax_, angMax_);
103  h_mResZMuMu_ = ZResHisto.make<TH1F>("ZToMuMuRecoMassResolution",
104  "Z Reco vs matched final state #mu #mu mass Difference (GeV/c^{2})",
105  nbinsMassRes_, -massResMax_, massResMax_);
106  h_mRatioZMuMu_ = ZResHisto.make<TH1F>("ZToMuMuRecoMassRatio",
107  "Z Reco vs matched final state #mu #mu mass Ratio",
108  4000, 0, 2);
109  h_mResZMuMuMC_ = ZResHisto.make<TH1F>("ZToMuMuMCMassResolution",
110  "Z vs final state #mu #mu MC mass Difference (GeV/c^{2})",
111  nbinsMassRes_/2 + 1, -2*massResMax_/nbinsMassRes_, massResMax_);
112  h_mRatioZMuMuMC_ = ZResHisto.make<TH1F>("ZToMuMuMCMassRatio",
113  "Z vs final state #mu #mu MC mass Ratio",
114  2002, 0.999, 2);
115 }
116 
117 void ZHistogrammer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
118  cout << ">>> Z Histogrammer analyze" << endl;
119  Handle<CandidateCollection> z;
120  Handle<CandidateCollection> gen;
121  Handle<CandMatchMap> match;
122  event.getByLabel(z_, z);
123  event.getByLabel(gen_, gen);
124  event.getByLabel(match_, match);
125  h_nZ_->Fill(z->size());
126  for(unsigned int i = 0; i < z->size(); ++i) {
127  const Candidate &zCand = (*z)[i];
128  h_mZ_->Fill(zCand.mass());
129  h_ptZ_->Fill(zCand.pt());
130  h_phiZ_->Fill(zCand.phi());
131  h_thetaZ_->Fill(zCand.theta());
132  h_etaZ_->Fill(zCand.eta());
133  h_rapidityZ_->Fill(zCand.rapidity());
134  CandidateRef zCandRef(z, i);
135  CandidateRef zMCMatch = (*match)[zCandRef];
136  if(zMCMatch.isNonnull() && zMCMatch->pdgId()==23) {
137  h_mResZ_->Fill(zCandRef->mass() - zMCMatch->mass());
138  h_ptResZ_->Fill(zCandRef->pt() - zMCMatch->pt());
139  h_phiResZ_->Fill(zCandRef->phi() - zMCMatch->phi());
140  h_thetaResZ_->Fill(zCandRef->theta() - zMCMatch->theta());
141  h_etaResZ_->Fill(zCandRef->eta() - zMCMatch->eta());
142  h_rapidityResZ_->Fill(zCandRef->rapidity() - zMCMatch->rapidity());
143  const Candidate * dau0 = zMCMatch->daughter(0);
144  const Candidate * dau1 = zMCMatch->daughter(1);
145  for(unsigned int i0 = 0; i0 < dau0->numberOfDaughters(); ++i0) {
146  const Candidate * ddau0 = dau0->daughter(i0);
147  if(abs(ddau0->pdgId())==13 && ddau0->status()==1) {
148  dau0 = ddau0; break;
149  }
150  }
151  for(unsigned int i1 = 0; i1 < dau1->numberOfDaughters(); ++i1) {
152  const Candidate * ddau1 = dau1->daughter(i1);
153  if(abs(ddau1->pdgId())==13 && ddau1->status()==1) {
154  dau1 = ddau1; break;
155  }
156  }
157  assert(abs(dau0->pdgId())==13 && dau0->status()==1);
158  assert(abs(dau1->pdgId())==13 && dau1->status()==1);
159  double invMass = (dau0->p4()+dau1->p4()).mass();
160  h_invmMuMu_->Fill(invMass);
161  h_mResZMuMu_->Fill(zCand.mass() - invMass);
162  h_mRatioZMuMu_->Fill(zCand.mass()/invMass);
163  }
164  }
165  h_nZMC_->Fill(gen->size());
166  for(unsigned int i = 0; i < gen->size(); ++i) {
167  const Candidate &genCand = (*gen)[i];
168  if((genCand.pdgId() == 23) && (genCand.status() == 2)) //this is an intermediate Z0
169  cout << ">>> intermediate Z0 found, with " << genCand.numberOfDaughters()
170  << " daughters" << endl;
171  if((genCand.pdgId() == 23)&&(genCand.status() == 3)) { //this is a Z0
172  cout << ">>> Z0 found, with " << genCand.numberOfDaughters()
173  << " daughters" << endl;
174  h_mZMC_->Fill(genCand.mass());
175  h_ptZMC_->Fill(genCand.pt());
176  h_phiZMC_->Fill(genCand.phi());
177  h_thetaZMC_->Fill(genCand.theta());
178  h_etaZMC_->Fill(genCand.eta());
179  h_rapidityZMC_->Fill(genCand.rapidity());
180  Particle::LorentzVector pZ(0, 0, 0, 0);
181  int nMu = 0;
182  for(unsigned int j = 0; j < genCand.numberOfDaughters(); ++j) {
183  const Candidate *dauGen = genCand.daughter(j);
184  /*
185  if((dauGen->pdgId() == 23) && (dauGen->status() == 2)) {
186  h_mZ2vs3MC_->Fill(genCand.mass() - dauGen->mass());
187  h_ptZ2vs3MC_->Fill(genCand.pt() - dauGen->pt());
188  h_phiZ2vs3MC_->Fill(genCand.phi() - dauGen->phi());
189  h_thetaZ2vs3MC_->Fill(genCand.theta() - dauGen->theta());
190  h_etaZ2vs3MC_->Fill(genCand.eta() - dauGen->eta());
191  h_rapidityZ2vs3MC_->Fill(genCand.rapidity() - dauGen->rapidity());
192  }
193  */
194  if((abs(dauGen->pdgId()) == 13) && (dauGen->numberOfDaughters() != 0)) {
195  //we are looking for photons of final state radiation
196  cout << ">>> The muon " << j
197  << " has " << dauGen->numberOfDaughters() << " daughters" <<endl;
198  for(unsigned int k = 0; k < dauGen->numberOfDaughters(); ++k) {
199  const Candidate * dauMuGen = dauGen->daughter(k);
200  cout << ">>> Mu " << j
201  << " daughter MC " << k
202  << " PDG Id " << dauMuGen->pdgId()
203  << ", status " << dauMuGen->status()
204  << ", charge " << dauMuGen->charge()
205  << endl;
206  if(abs(dauMuGen->pdgId()) == 13 && dauMuGen->status() ==1) {
207  pZ += dauMuGen->p4();
208  nMu ++;
209  }
210  }
211  }
212  }
213  assert(nMu == 2);
214  double mZ = pZ.mass();
215  h_invmMuMuMC_->Fill(mZ);
216  h_mResZMuMuMC_->Fill(genCand.mass() - mZ);
217  h_mRatioZMuMuMC_->Fill(genCand.mass()/mZ);
218  }
219  }
220 }
221 
222 #include "FWCore/Framework/interface/MakerMacros.h"
223 
224 DEFINE_FWK_MODULE(ZHistogrammer);
225 
i
int i
Definition: DBlmapReader.cc:9
reco::Candidate::daughter
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
ZHistogrammer::h_mResZMuMuMC_
TH1F * h_mResZMuMuMC_
Definition: ZHistogrammer.cc:20
ZHistogrammer::h_mRatioZMuMu_
TH1F * h_mRatioZMuMu_
Definition: ZHistogrammer.cc:19
reco::Candidate::mass
virtual float mass() const =0
mass
ZHistogrammer::nbinsPt_
unsigned int nbinsPt_
Definition: ZHistogrammer.cc:11
reco::Candidate::eta
virtual float eta() const =0
momentum pseudorapidity
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
reco::Candidate::status
virtual int status() const =0
status word
abs
#define abs(x)
Definition: mlp_lapack.h:159
reco::Candidate::phi
virtual float phi() const =0
momentum azimuthal angle
reco::Candidate::rapidity
virtual double rapidity() const =0
rapidity
ZHistogrammer::nbinsAng_
unsigned int nbinsAng_
Definition: ZHistogrammer.cc:11
detailsBasic3DVector::z
float float float z
Definition: extBasic3DVector.h:15
ZHistogrammer::h_rapidityZ_
TH1F * h_rapidityZ_
Definition: ZHistogrammer.cc:13
relval_steps.gen
def gen
Definition: relval_steps.py:239
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
reco::Candidate::numberOfDaughters
virtual size_type numberOfDaughters() const =0
number of daughters
reco::Candidate::pt
virtual float pt() const =0
transverse momentum
reco::Candidate::theta
virtual double theta() const =0
momentum polar angle
ZHistogrammer::h_invmMuMuMC_
TH1F * h_invmMuMuMC_
Definition: ZHistogrammer.cc:16
L1Trigger_dataformats.reco
dictionary reco
Definition: L1Trigger_dataformats.py:9
ZHistogrammer::analyze
virtual void analyze(const edm::Event &event, const edm::EventSetup &setup)
Definition: ZHistogrammer.cc:117
ZHistogrammer::h_rapidityZMC_
TH1F * h_rapidityZMC_
Definition: ZHistogrammer.cc:15
ZHistogrammer::h_thetaZMC_
TH1F * h_thetaZMC_
Definition: ZHistogrammer.cc:15
j
int j
Definition: DBlmapReader.cc:9
ZHistogrammer::ZHistogrammer
ZHistogrammer(const edm::ParameterSet &pset)
Definition: ZHistogrammer.cc:40
ZHistogrammer::h_rapidityResZ_
TH1F * h_rapidityResZ_
Definition: ZHistogrammer.cc:18
reco::Candidate::charge
virtual int charge() const =0
electric charge
event
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
ZHistogrammer::nbinsMassRes_
unsigned int nbinsMassRes_
Definition: ZHistogrammer.cc:11
ZHistogrammer::z_
edm::InputTag z_
Definition: ZHistogrammer.cc:10
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
reco::Candidate::pdgId
virtual int pdgId() const =0
PDG identifier.
TFileDirectory::mkdir
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileDirectory.cc:162
ZHistogrammer::match_
edm::InputTag match_
Definition: ZHistogrammer.cc:10
ZHistogrammer::h_mResZMuMu_
TH1F * h_mResZMuMu_
Definition: ZHistogrammer.cc:19
ZHistogrammer::h_invmMuMu_
TH1F * h_invmMuMu_
Definition: ZHistogrammer.cc:14
TFileDirectory::make
T * make() const
make new ROOT object
Definition: TFileDirectory.h:55
ZHistogrammer::massResMax_
double massResMax_
Definition: ZHistogrammer.cc:12
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
match
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
ZHistogrammer::h_mRatioZMuMuMC_
TH1F * h_mRatioZMuMuMC_
Definition: ZHistogrammer.cc:20
ZHistogrammer::nbinsMass_
unsigned int nbinsMass_
Definition: ZHistogrammer.cc:11
reco::Particle::LorentzVector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:25
HcalObjRepresent::setup
void setup(std::vector< TH2F > &depth, std::string name, std::string units="")
Definition: HcalObjRepresent.h:628
ZHistogrammer::gen_
edm::InputTag gen_
Definition: ZHistogrammer.cc:10
reco::Candidate::p4
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
ZHistogrammer::h_thetaResZ_
TH1F * h_thetaResZ_
Definition: ZHistogrammer.cc:18
Generated for CMSSW Reference Manual by doxygen 1.8.5