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.
6 #include "TH1.h"
7 
8 class ZHistogrammer : public edm::EDAnalyzer {
9 public:
11 private:
12  void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
19  TH1F *h_invmMuMu_;
22  //TH1F *h_mZ2vs3MC_, *h_ptZ2vs3MC_, *h_phiZ2vs3MC_, *h_thetaZ2vs3MC_, *h_etaZ2vs3MC_, *h_rapidityZ2vs3MC_;
26 };
27 
35 #include <cmath>
36 #include <iostream>
37 
38 using namespace std;
39 using namespace reco;
40 using namespace edm;
41 
43  zToken_(consumes<CandidateCollection>(pset.getParameter<InputTag>("z"))),
44  genToken_(consumes<CandidateCollection>(pset.getParameter<InputTag>("gen"))),
45  matchToken_(consumes<CandMatchMap>(pset.getParameter<InputTag>("match"))),
46  nbinsMass_(pset.getUntrackedParameter<unsigned int>("nbinsMass")),
47  nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")),
48  nbinsAng_(pset.getUntrackedParameter<unsigned int>("nbinsAng")),
49  nbinsMassRes_(pset.getUntrackedParameter<unsigned int>("nbinsMassRes")),
50  massMax_(pset.getUntrackedParameter<double>("massMax")),
51  ptMax_(pset.getUntrackedParameter<double>("ptMax")),
52  angMax_(pset.getUntrackedParameter<double>("angMax")),
53  massResMax_(pset.getUntrackedParameter<double>("massResMax")) {
54  cout << ">>> Z Histogrammer constructor" << endl;
56  TFileDirectory ZHisto = fs->mkdir( "ZRecoHisto" );
57  TFileDirectory ZMCHisto = fs->mkdir( "ZMCHisto" );
58  TFileDirectory ZResHisto = fs->mkdir( "ZResolutionHisto" );
59  //TFileDirectory Z2vs3MCHisto = fs->mkdir( "Z2vs3MCHisto" );
60  h_nZ_ = ZHisto.make<TH1F>("ZNumber", "number of Z particles", 11, -0.5, 10.5);
61  h_mZ_ = ZHisto.make<TH1F>("ZMass", "Z mass (GeV/c^{2})", nbinsMass_, 0, massMax_);
62  h_ptZ_ = ZHisto.make<TH1F>("ZPt", "Z p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
63  h_phiZ_ = ZHisto.make<TH1F>("ZPhi", "Z #phi", nbinsAng_, -angMax_, angMax_);
64  h_thetaZ_ = ZHisto.make<TH1F>("Ztheta", "Z #theta", nbinsAng_, 0, angMax_);
65  h_etaZ_ = ZHisto.make<TH1F>("ZEta", "Z #eta", nbinsAng_, -angMax_, angMax_);
66  h_rapidityZ_ = ZHisto.make<TH1F>("ZRapidity", "Z rapidity", nbinsAng_, -angMax_, angMax_);
67  h_invmMuMu_ = ZHisto.make<TH1F>("MuMuMass", "#mu #mu invariant mass",
68  nbinsMass_, 0, massMax_);
69  h_nZMC_ = ZMCHisto.make<TH1F>("ZMCNumber", "number of Z MC particles", 11, -0.5, 10.5);
70  h_mZMC_ = ZMCHisto.make<TH1F>("ZMCMass", "Z MC mass (GeV/c^{2})", nbinsMass_, 0, massMax_);
71  h_ptZMC_ = ZMCHisto.make<TH1F>("ZMCPt", "Z MC p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
72  h_phiZMC_ = ZMCHisto.make<TH1F>("ZMCPhi", "Z MC #phi", nbinsAng_, -angMax_, angMax_);
73  h_thetaZMC_ = ZMCHisto.make<TH1F>("ZMCTheta", "Z MC #theta", nbinsAng_, 0, angMax_);
74  h_etaZMC_ = ZMCHisto.make<TH1F>("ZMCEta", "Z MC #eta", nbinsAng_, -angMax_, angMax_);
75  h_rapidityZMC_ = ZMCHisto.make<TH1F>("ZMCRapidity", "Z MC rapidity",
77  h_invmMuMuMC_ = ZMCHisto.make<TH1F>("MuMuMCMass", "#mu #mu MC invariant mass",
78  nbinsMass_, 0, massMax_);
79  /*
80  h_mZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCMass", "Z MC st 2 vs st 3 mass (GeV/c^{2})",
81  nbinsMassRes_, -massResMax_, massResMax_);
82  h_ptZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCPt", "Z MC st 2 vs st 3 p_{t} (GeV/c)",
83  nbinsPt_, -ptMax_, ptMax_);
84  h_phiZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCPhi", "Z MC st 2 vs st 3 #phi",
85  nbinsAng_, -angMax_, angMax_);
86  h_thetaZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCTheta", "Z MC st 2 vs st 3 #theta",
87  nbinsAng_, -angMax_, angMax_);
88  h_etaZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCEta", "Z MC st 2 vs st 3 #eta",
89  nbinsAng_, -angMax_, angMax_);
90  h_rapidityZ2vs3MC_ = Z2vs3MCHisto.make<TH1F>("Z2vs3MCRapidity", "Z MC st 2 vs st 3 rapidity",
91  nbinsAng_, -angMax_, angMax_);
92  */
93  h_mResZ_ = ZResHisto.make<TH1F>("ZMassResolution", "Z mass Resolution (GeV/c^{2})",
95  h_ptResZ_ = ZResHisto.make<TH1F>("ZPtResolution", "Z p_{t} Resolution (GeV/c)",
97  h_phiResZ_ = ZResHisto.make<TH1F>("ZPhiResolution", "Z #phi Resolution",
99  h_thetaResZ_ = ZResHisto.make<TH1F>("ZThetaResolution", "Z #theta Resolution",
101  h_etaResZ_ = ZResHisto.make<TH1F>("ZEtaResolution", "Z #eta Resolution",
103  h_rapidityResZ_ = ZResHisto.make<TH1F>("ZRapidityResolution", "Z rapidity Resolution",
105  h_mResZMuMu_ = ZResHisto.make<TH1F>("ZToMuMuRecoMassResolution",
106  "Z Reco vs matched final state #mu #mu mass Difference (GeV/c^{2})",
108  h_mRatioZMuMu_ = ZResHisto.make<TH1F>("ZToMuMuRecoMassRatio",
109  "Z Reco vs matched final state #mu #mu mass Ratio",
110  4000, 0, 2);
111  h_mResZMuMuMC_ = ZResHisto.make<TH1F>("ZToMuMuMCMassResolution",
112  "Z vs final state #mu #mu MC mass Difference (GeV/c^{2})",
113  nbinsMassRes_/2 + 1, -2*massResMax_/nbinsMassRes_, massResMax_);
114  h_mRatioZMuMuMC_ = ZResHisto.make<TH1F>("ZToMuMuMCMassRatio",
115  "Z vs final state #mu #mu MC mass Ratio",
116  2002, 0.999, 2);
117 }
118 
120  cout << ">>> Z Histogrammer analyze" << endl;
124  event.getByToken(zToken_, z);
125  event.getByToken(genToken_, gen);
126  event.getByToken(matchToken_, match);
127  h_nZ_->Fill(z->size());
128  for(unsigned int i = 0; i < z->size(); ++i) {
129  const Candidate &zCand = (*z)[i];
130  h_mZ_->Fill(zCand.mass());
131  h_ptZ_->Fill(zCand.pt());
132  h_phiZ_->Fill(zCand.phi());
133  h_thetaZ_->Fill(zCand.theta());
134  h_etaZ_->Fill(zCand.eta());
135  h_rapidityZ_->Fill(zCand.rapidity());
136  CandidateRef zCandRef(z, i);
137  CandidateRef zMCMatch = (*match)[zCandRef];
138  if(zMCMatch.isNonnull() && zMCMatch->pdgId()==23) {
139  h_mResZ_->Fill(zCandRef->mass() - zMCMatch->mass());
140  h_ptResZ_->Fill(zCandRef->pt() - zMCMatch->pt());
141  h_phiResZ_->Fill(zCandRef->phi() - zMCMatch->phi());
142  h_thetaResZ_->Fill(zCandRef->theta() - zMCMatch->theta());
143  h_etaResZ_->Fill(zCandRef->eta() - zMCMatch->eta());
144  h_rapidityResZ_->Fill(zCandRef->rapidity() - zMCMatch->rapidity());
145  const Candidate * dau0 = zMCMatch->daughter(0);
146  const Candidate * dau1 = zMCMatch->daughter(1);
147  for(unsigned int i0 = 0; i0 < dau0->numberOfDaughters(); ++i0) {
148  const Candidate * ddau0 = dau0->daughter(i0);
149  if(abs(ddau0->pdgId())==13 && ddau0->status()==1) {
150  dau0 = ddau0; break;
151  }
152  }
153  for(unsigned int i1 = 0; i1 < dau1->numberOfDaughters(); ++i1) {
154  const Candidate * ddau1 = dau1->daughter(i1);
155  if(abs(ddau1->pdgId())==13 && ddau1->status()==1) {
156  dau1 = ddau1; break;
157  }
158  }
159  assert(abs(dau0->pdgId())==13 && dau0->status()==1);
160  assert(abs(dau1->pdgId())==13 && dau1->status()==1);
161  double invMass = (dau0->p4()+dau1->p4()).mass();
162  h_invmMuMu_->Fill(invMass);
163  h_mResZMuMu_->Fill(zCand.mass() - invMass);
164  h_mRatioZMuMu_->Fill(zCand.mass()/invMass);
165  }
166  }
167  h_nZMC_->Fill(gen->size());
168  for(unsigned int i = 0; i < gen->size(); ++i) {
169  const Candidate &genCand = (*gen)[i];
170  if((genCand.pdgId() == 23) && (genCand.status() == 2)) //this is an intermediate Z0
171  cout << ">>> intermediate Z0 found, with " << genCand.numberOfDaughters()
172  << " daughters" << endl;
173  if((genCand.pdgId() == 23)&&(genCand.status() == 3)) { //this is a Z0
174  cout << ">>> Z0 found, with " << genCand.numberOfDaughters()
175  << " daughters" << endl;
176  h_mZMC_->Fill(genCand.mass());
177  h_ptZMC_->Fill(genCand.pt());
178  h_phiZMC_->Fill(genCand.phi());
179  h_thetaZMC_->Fill(genCand.theta());
180  h_etaZMC_->Fill(genCand.eta());
181  h_rapidityZMC_->Fill(genCand.rapidity());
182  Particle::LorentzVector pZ(0, 0, 0, 0);
183  int nMu = 0;
184  for(unsigned int j = 0; j < genCand.numberOfDaughters(); ++j) {
185  const Candidate *dauGen = genCand.daughter(j);
186  /*
187  if((dauGen->pdgId() == 23) && (dauGen->status() == 2)) {
188  h_mZ2vs3MC_->Fill(genCand.mass() - dauGen->mass());
189  h_ptZ2vs3MC_->Fill(genCand.pt() - dauGen->pt());
190  h_phiZ2vs3MC_->Fill(genCand.phi() - dauGen->phi());
191  h_thetaZ2vs3MC_->Fill(genCand.theta() - dauGen->theta());
192  h_etaZ2vs3MC_->Fill(genCand.eta() - dauGen->eta());
193  h_rapidityZ2vs3MC_->Fill(genCand.rapidity() - dauGen->rapidity());
194  }
195  */
196  if((abs(dauGen->pdgId()) == 13) && (dauGen->numberOfDaughters() != 0)) {
197  //we are looking for photons of final state radiation
198  cout << ">>> The muon " << j
199  << " has " << dauGen->numberOfDaughters() << " daughters" <<endl;
200  for(unsigned int k = 0; k < dauGen->numberOfDaughters(); ++k) {
201  const Candidate * dauMuGen = dauGen->daughter(k);
202  cout << ">>> Mu " << j
203  << " daughter MC " << k
204  << " PDG Id " << dauMuGen->pdgId()
205  << ", status " << dauMuGen->status()
206  << ", charge " << dauMuGen->charge()
207  << endl;
208  if(abs(dauMuGen->pdgId()) == 13 && dauMuGen->status() ==1) {
209  pZ += dauMuGen->p4();
210  nMu ++;
211  }
212  }
213  }
214  }
215  assert(nMu == 2);
216  double mZ = pZ.mass();
217  h_invmMuMuMC_->Fill(mZ);
218  h_mResZMuMuMC_->Fill(genCand.mass() - mZ);
219  h_mRatioZMuMuMC_->Fill(genCand.mass()/mZ);
220  }
221  }
222 }
223 
225 
227 
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:251
TH1F * h_mResZMuMuMC_
TH1F * h_mRatioZMuMu_
void analyze(const edm::Event &event, const edm::EventSetup &setup) override
unsigned int nbinsPt_
virtual const Candidate * daughter(size_type i) const =0
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
size_type size() const
Definition: OwnVector.h:264
def setup(process, global_tag, zero_tesla=False)
Definition: GeneralSetup.py:2
unsigned int nbinsAng_
TH1F * h_rapidityZ_
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
virtual int status() const =0
status word
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
virtual const LorentzVector & p4() const =0
four-momentum Lorentz vector
TH1F * h_invmMuMuMC_
virtual int pdgId() const =0
PDG identifier.
TH1F * h_rapidityZMC_
TH1F * h_thetaZMC_
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ZHistogrammer(const edm::ParameterSet &pset)
TH1F * h_rapidityResZ_
T * make(const Args &...args) const
make new ROOT object
def gen(fragment, howMuch)
Production test section ####.
unsigned int nbinsMassRes_
int k[5][pyjets_maxn]
virtual double theta() const =0
momentum polar angle
virtual double eta() const =0
momentum pseudorapidity
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
virtual double pt() const =0
transverse momentum
edm::EDGetTokenT< reco::CandidateCollection > zToken_
virtual double mass() const =0
mass
edm::EDGetTokenT< reco::CandidateCollection > genToken_
TH1F * h_mResZMuMu_
TH1F * h_invmMuMu_
virtual double rapidity() const =0
rapidity
virtual int charge() const =0
electric charge
fixed size matrix
HLT enums.
edm::EDGetTokenT< reco::CandMatchMap > matchToken_
double massResMax_
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:10
TH1F * h_mRatioZMuMuMC_
unsigned int nbinsMass_
virtual size_type numberOfDaughters() const =0
number of daughters
virtual double phi() const =0
momentum azimuthal angle
Definition: event.py:1
math::PtEtaPhiELorentzVectorF LorentzVector
TH1F * h_thetaResZ_