CMS 3D CMS Logo
zPdfUnc.cc
Go to the documentation of this file.
1 #include "FWCore/Framework/interface/EDAnalyzer.h"
2 #include "FWCore/Utilities/interface/InputTag.h"
3 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
4 #include <vector>
5 #include "TH1.h"
6 
7 class zPdfUnc : public edm::EDAnalyzer {
8 public:
9  zPdfUnc(const edm::ParameterSet& pset);
10 
11 private:
12  void analyze(const edm::Event& event, const edm::EventSetup& setup) override;
13  void endJob() override;
14  edm::EDGetTokenT<reco::GenParticleCollection> genToken_;
15  edm::EDGetTokenT<std::vector<double> > pdfweightsToken_;
16  unsigned int pdfmember_, nbinsMass_, nbinsPt_, nbinsAng_;
17  double massMax_, ptMax_, angMax_;
18  double accPtMin_, accMassMin_, accMassMax_, accEtaMin_, accEtaMax_, accMassMinDenominator_;
19  TH1F* h_nZ_;
20  TH1F *h_mZMC_, *h_ptZMC_, *h_phiZMC_, *h_thetaZMC_, *h_etaZMC_, *h_rapidityZMC_;
21  TH1F *hardpt, *softpt, *hardeta, *softeta;
22  TH1F* h_weight_histo;
23  bool isMCatNLO_;
24  double nAcc_, nAccReW_, nBothMuHasZHasGrandMa_;
25  std::string filename_;
26 };
27 
28 #include "DataFormats/Candidate/interface/Candidate.h"
29 #include "DataFormats/Candidate/interface/CandidateFwd.h"
30 #include "DataFormats/Candidate/interface/CandMatchMap.h"
31 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
32 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
33 
34 #include "HepMC/WeightContainer.h"
35 #include "HepMC/GenEvent.h"
36 #include "FWCore/Framework/interface/Event.h"
37 #include "DataFormats/Common/interface/Handle.h"
38 #include "FWCore/ParameterSet/interface/ParameterSet.h"
39 #include "FWCore/ServiceRegistry/interface/Service.h"
40 #include "CommonTools/UtilAlgos/interface/TFileService.h"
41 #include <cmath>
42 #include <iostream>
43 #include <fstream>
44 
45 using namespace std;
46 using namespace reco;
47 using namespace edm;
48 
49 zPdfUnc::zPdfUnc(const ParameterSet& pset)
50  : genToken_(consumes<GenParticleCollection>(pset.getParameter<InputTag>("genParticles"))),
51  pdfweightsToken_(consumes<std::vector<double> >(pset.getParameter<InputTag>("pdfweights"))),
52  pdfmember_(pset.getUntrackedParameter<unsigned int>("pdfmember")),
53  nbinsMass_(pset.getUntrackedParameter<unsigned int>("nbinsMass")),
54  nbinsPt_(pset.getUntrackedParameter<unsigned int>("nbinsPt")),
55  nbinsAng_(pset.getUntrackedParameter<unsigned int>("nbinsAng")),
56  massMax_(pset.getUntrackedParameter<double>("massMax")),
57  ptMax_(pset.getUntrackedParameter<double>("ptMax")),
58  angMax_(pset.getUntrackedParameter<double>("angMax")),
59  accPtMin_(pset.getUntrackedParameter<double>("accPtMin")),
60  accMassMin_(pset.getUntrackedParameter<double>("accMassMin")),
61  accMassMax_(pset.getUntrackedParameter<double>("accMassMax")),
62  accEtaMin_(pset.getUntrackedParameter<double>("accEtaMin")),
63  accEtaMax_(pset.getUntrackedParameter<double>("accEtaMax")),
64  accMassMinDenominator_(pset.getUntrackedParameter<double>("accMassMinDenominator")),
65  isMCatNLO_(pset.getUntrackedParameter<bool>("isMCatNLO")),
66  filename_(pset.getUntrackedParameter<std::string>("outfilename")) {
67  cout << ">>> Z Histogrammer constructor" << endl;
68  Service<TFileService> fs;
69 
70  TFileDirectory ZMCHisto = fs->mkdir("ZMCHisto");
71  h_nZ_ = ZMCHisto.make<TH1F>("ZNumber", "number of Z particles", 11, -0.5, 10.5);
72  h_weight_histo = ZMCHisto.make<TH1F>("weight_histo", "weight_histo", 20, -10, 10);
73 
74  h_mZMC_ = ZMCHisto.make<TH1F>("ZMCMass", "Z MC mass (GeV/c^{2})", nbinsMass_, 0, massMax_);
75  h_ptZMC_ = ZMCHisto.make<TH1F>("ZMCPt", "Z MC p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
76  hardpt = ZMCHisto.make<TH1F>("hardpt", "hard muon p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
77  softpt = ZMCHisto.make<TH1F>("softpt", "soft muon p_{t} (GeV/c)", nbinsPt_, 0, ptMax_);
78 
79  h_phiZMC_ = ZMCHisto.make<TH1F>("ZMCPhi", "Z MC #phi", nbinsAng_, -angMax_, angMax_);
80  h_thetaZMC_ = ZMCHisto.make<TH1F>("ZMCTheta", "Z MC #theta", nbinsAng_, 0, angMax_);
81  h_etaZMC_ = ZMCHisto.make<TH1F>("ZMCEta", "Z MC #eta", nbinsAng_, -angMax_, angMax_);
82  h_rapidityZMC_ = ZMCHisto.make<TH1F>("ZMCRapidity", "Z MC y", nbinsAng_, -angMax_, angMax_);
83 
84  hardeta = ZMCHisto.make<TH1F>("hard muon eta", "hard muon #eta", nbinsAng_, -angMax_, angMax_);
85  softeta = ZMCHisto.make<TH1F>("soft muon eta", "soft muon #eta", nbinsAng_, -angMax_, angMax_);
86  nAcc_ = 0.;
87  nAccReW_ = 0;
88  nBothMuHasZHasGrandMa_ = 0;
89 }
90 
91 void zPdfUnc::analyze(const edm::Event& event, const edm::EventSetup& setup) {
92  cout << ">>> Z Histogrammer analyze" << endl;
93 
94  Handle<GenParticleCollection> gen;
95  Handle<vector<double> > pdfweights;
96 
97  event.getByToken(genToken_, gen);
98  event.getByToken(pdfweightsToken_, pdfweights);
99 
100  // get weight and fill it to histogram
101  std::vector<double> weights = (*pdfweights);
102 
103  unsigned int nmembers = weights.size();
104 
105  double weight = weights[pdfmember_];
106 
107  h_weight_histo->Fill(weight);
108 
109  cout << "found nmember: " << nmembers << endl;
110  cout << "taken the member n " << pdfmember_ << " weight= " << weight << endl;
111 
112  std::vector<GenParticle> muons;
113 
114  for (unsigned int i = 0; i < gen->size(); ++i) {
115  const GenParticle& muMC = (*gen)[i];
116  // filling only muons coming form Z
117  if (abs(muMC.pdgId()) == 13 && muMC.status() == 1 && muMC.numberOfMothers() > 0) {
118  if (muMC.mother()->numberOfMothers() > 0) {
119  cout << "I'm getting a muon \n"
120  << "with "
121  << "muMC.numberOfMothers() " << muMC.numberOfMothers() << "\n the first mother has pdgId "
122  << muMC.mother()->pdgId() << "with "
123  << "muMC.mother()->numberOfMothers() " << muMC.mother()->numberOfMothers()
124  << "\n the first grandma has pdgId " << muMC.mother()->mother()->pdgId() << endl;
125  if (muMC.mother()->mother()->pdgId() == 23)
126  muons.push_back(muMC);
127  }
128  }
129  }
130 
131  cout << "finally I selected " << muons.size() << " muons" << endl;
132 
133  // if there are at least two muons,
134  // calculate invarant mass of first two and fill it into histogram
135 
136  double inv_mass = 0.0;
137  double Zpt_ = 0.0;
138  double Zeta_ = 0.0;
139  double Ztheta_ = 0.0;
140  double Zphi_ = 0.0;
141  double Zrapidity_ = 0.0;
142 
143  if (muons.size() > 1) {
144  if (muons[0].mother()->mother()->pdgId() == 23 && muons[1].mother()->mother()->pdgId() == 23)
145  nBothMuHasZHasGrandMa_++;
146  math::XYZTLorentzVector tot_momentum(muons[0].p4());
147  math::XYZTLorentzVector mom2(muons[1].p4());
148  tot_momentum += mom2;
149  inv_mass = sqrt(tot_momentum.mass2());
150 
151  Zpt_ = tot_momentum.pt();
152  Zeta_ = tot_momentum.eta();
153  Ztheta_ = tot_momentum.theta();
154  Zphi_ = tot_momentum.phi();
155  Zrapidity_ = tot_momentum.Rapidity();
156 
157  // IMPORTANT: use the weight of the event ...
158 
159  double weight_sign = weight;
160  //double weight_sign = 1. ;
161  // fill the denominator numbers only if mass>accMassMinDenominator and go on in that case
162  if (inv_mass > accMassMinDenominator_) {
163  h_mZMC_->Fill(inv_mass, weight_sign);
164  h_ptZMC_->Fill(Zpt_, weight_sign);
165  h_etaZMC_->Fill(Zeta_, weight_sign);
166  h_thetaZMC_->Fill(Ztheta_, weight_sign);
167  h_phiZMC_->Fill(Zphi_, weight_sign);
168  h_rapidityZMC_->Fill(Zrapidity_, weight_sign);
169 
170  double pt1 = muons[0].pt();
171  double pt2 = muons[1].pt();
172  double eta1 = muons[0].eta();
173  double eta2 = muons[1].eta();
174  if (pt1 > pt2) {
175  hardpt->Fill(pt1, weight_sign);
176  softpt->Fill(pt2, weight_sign);
177  hardeta->Fill(eta1, weight_sign);
178  softeta->Fill(eta2, weight_sign);
179  } else {
180  hardpt->Fill(pt2, weight_sign);
181  softpt->Fill(pt1, weight_sign);
182  hardeta->Fill(eta2, weight_sign);
183  softeta->Fill(eta1, weight_sign);
184  }
185 
186  //evaluating the geometric acceptance
187  if (pt1 >= accPtMin_ && pt2 >= accPtMin_ && fabs(eta1) >= accEtaMin_ && fabs(eta2) >= accEtaMin_ &&
188  fabs(eta1) <= accEtaMax_ && fabs(eta2) <= accEtaMax_ && inv_mass >= accMassMin_ && inv_mass <= accMassMax_) {
189  nAcc_++;
190  nAccReW_ += weight;
191  }
192  }
193  }
194 }
195 
196 void zPdfUnc::endJob() {
197  cout << " number of events accepted :" << nAcc_ << endl;
198  cout << " number of events accepted reweigthed :" << nAccReW_ << endl;
199  double nev = h_mZMC_->GetEntries();
200  double nev_weigthed = h_mZMC_->Integral();
201  cout << " number of total events :" << nev << endl;
202  cout << " number of total weighted events :" << nev_weigthed << endl;
203  cout << " number of cases in which BothMuHasZHasGrandMa :" << nBothMuHasZHasGrandMa_ << endl;
204  double eff = (double)nAcc_ / (double)h_mZMC_->GetEntries();
205  double eff_rew = (double)nAccReW_ / (double)h_mZMC_->Integral(0, nbinsMass_ + 1);
206  double err = sqrt(eff * (1. - eff) / (double)h_mZMC_->GetEntries());
207  double err_rew = sqrt(eff_rew * (1. - eff_rew) / (double)h_mZMC_->Integral(0, nbinsMass_ + 1));
208  cout << " geometric acceptance: " << eff << "+/-" << err << endl;
209  cout << " geometric acceptance reweighted: " << eff_rew << "+/-" << err_rew << endl;
210 
211  ofstream myfile;
212  myfile.open(filename_.c_str(), std::ios::app);
213  myfile << eff << " " << eff_rew << " " << nev << " " << nev_weigthed << endl;
214  myfile.close();
215 }
216 #include "FWCore/Framework/interface/MakerMacros.h"
217 
218 DEFINE_FWK_MODULE(zPdfUnc);
PDWG_BPHSkim_cff.muons
muons
Definition: PDWG_BPHSkim_cff.py:47
zPdfUnc::endJob
void endJob() override
Definition: zPdfUnc.cc:196
Handle.h
zPdfUnc::nbinsAng_
unsigned int nbinsAng_
Definition: zPdfUnc.cc:16
zPdfUnc::h_ptZMC_
TH1F * h_ptZMC_
Definition: zPdfUnc.cc:20
electrons_cff.bool
bool
Definition: electrons_cff.py:393
mps_fire.i
i
Definition: mps_fire.py:428
reco::LeafCandidate::status
int status() const final
status word
Definition: LeafCandidate.h:180
edm::EDGetTokenT< reco::GenParticleCollection >
reco::CompositeRefCandidateT::mother
const Candidate * mother(size_type=0) const override
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode)
TFileDirectory::make
T * make(const Args &... args) const
make new ROOT object
Definition: TFileDirectory.h:53
edm
HLT enums.
Definition: AlignableModifier.h:19
mps_merge.weight
weight
Definition: mps_merge.py:88
gather_cfg.cout
cout
Definition: gather_cfg.py:144
reco::GenParticleCollection
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
Definition: GenParticleFwd.h:13
zPdfUnc::isMCatNLO_
bool isMCatNLO_
Definition: zPdfUnc.cc:23
reco::CompositeRefCandidateT::numberOfMothers
size_t numberOfMothers() const override
number of mothers
zPdfUnc
Definition: zPdfUnc.cc:7
EDAnalyzer.h
zPdfUnc::h_rapidityZMC_
TH1F * h_rapidityZMC_
Definition: zPdfUnc.cc:20
TFileDirectory
Definition: TFileDirectory.h:24
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
reco::Candidate::mother
virtual const Candidate * mother(size_type i=0) const =0
return pointer to mother
edm::Handle
Definition: AssociativeIterator.h:50
GenParticle
Definition: GenParticle.py:1
zPdfUnc::massMax_
double massMax_
Definition: zPdfUnc.cc:17
singleTopDQM_cfi.setup
setup
Definition: singleTopDQM_cfi.py:37
CandMatchMap.h
edm::EDAnalyzer
Definition: EDAnalyzer.h:28
zPdfUnc::pdfmember_
unsigned int pdfmember_
Definition: zPdfUnc.cc:16
zPdfUnc::accMassMinDenominator_
double accMassMinDenominator_
Definition: zPdfUnc.cc:18
GenParticle.h
CandidateFwd.h
MakerMacros.h
reco::Candidate::numberOfMothers
virtual size_type numberOfMothers() const =0
number of mothers (zero or one in most of but not all the cases)
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Service.h
HLT_FULL_cff.weights
weights
Definition: HLT_FULL_cff.py:99178
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
zPdfUnc::nAccReW_
double nAccReW_
Definition: zPdfUnc.cc:24
HLT_FULL_cff.eta2
eta2
Definition: HLT_FULL_cff.py:9542
GenParticleFwd.h
HLT_FULL_cff.pt1
pt1
Definition: HLT_FULL_cff.py:9870
gen
Definition: PythiaDecays.h:13
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
TFileService.h
HLT_FULL_cff.eta1
eta1
Definition: HLT_FULL_cff.py:9541
TFileService::mkdir
TFileDirectory mkdir(const std::string &dir, const std::string &descr="")
create a new subdirectory
Definition: TFileService.h:69
zPdfUnc::filename_
std::string filename_
Definition: zPdfUnc.cc:25
zPdfUnc::h_nZ_
TH1F * h_nZ_
Definition: zPdfUnc.cc:19
edm::ParameterSet
Definition: ParameterSet.h:47
zPdfUnc::softpt
TH1F * softpt
Definition: zPdfUnc.cc:21
Event.h
zPdfUnc::ptMax_
double ptMax_
Definition: zPdfUnc.cc:17
zPdfUnc::h_weight_histo
TH1F * h_weight_histo
Definition: zPdfUnc.cc:22
reco::LeafCandidate::pdgId
int pdgId() const final
PDG identifier.
Definition: LeafCandidate.h:176
zPdfUnc::hardpt
TH1F * hardpt
Definition: zPdfUnc.cc:21
edm::Service< TFileService >
createfilelist.int
int
Definition: createfilelist.py:10
p4
double p4[4]
Definition: TauolaWrapper.h:92
trackerHitRTTI::vector
Definition: trackerHitRTTI.h:21
zPdfUnc::nbinsMass_
unsigned int nbinsMass_
Definition: zPdfUnc.cc:16
zPdfUnc::analyze
void analyze(const edm::Event &event, const edm::EventSetup &setup) override
Definition: zPdfUnc.cc:91
zPdfUnc::accEtaMax_
double accEtaMax_
Definition: zPdfUnc.cc:18
edm::EventSetup
Definition: EventSetup.h:57
HLT_FULL_cff.pt2
pt2
Definition: HLT_FULL_cff.py:9872
zPdfUnc::pdfweightsToken_
edm::EDGetTokenT< std::vector< double > > pdfweightsToken_
Definition: zPdfUnc.cc:15
submitPVResolutionJobs.err
err
Definition: submitPVResolutionJobs.py:85
zPdfUnc::accMassMax_
double accMassMax_
Definition: zPdfUnc.cc:18
reco::Candidate::pdgId
virtual int pdgId() const =0
PDG identifier.
zPdfUnc::genToken_
edm::EDGetTokenT< reco::GenParticleCollection > genToken_
Definition: zPdfUnc.cc:14
InputTag.h
zPdfUnc::angMax_
double angMax_
Definition: zPdfUnc.cc:17
zPdfUnc::h_phiZMC_
TH1F * h_phiZMC_
Definition: zPdfUnc.cc:20
std
Definition: JetResolutionObject.h:76
zPdfUnc::nAcc_
double nAcc_
Definition: zPdfUnc.cc:24
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
zPdfUnc::h_etaZMC_
TH1F * h_etaZMC_
Definition: zPdfUnc.cc:20
zPdfUnc::nBothMuHasZHasGrandMa_
double nBothMuHasZHasGrandMa_
Definition: zPdfUnc.cc:24
relval_steps.gen
def gen(fragment, howMuch)
Production test section ####.
Definition: relval_steps.py:509
zPdfUnc::accEtaMin_
double accEtaMin_
Definition: zPdfUnc.cc:18
Candidate.h
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
ParameterSet.h
HepMCProduct.h
event
Definition: event.py:1
edm::Event
Definition: Event.h:73
zPdfUnc::accMassMin_
double accMassMin_
Definition: zPdfUnc.cc:18
zPdfUnc::nbinsPt_
unsigned int nbinsPt_
Definition: zPdfUnc.cc:16
zPdfUnc::softeta
TH1F * softeta
Definition: zPdfUnc.cc:21
zPdfUnc::hardeta
TH1F * hardeta
Definition: zPdfUnc.cc:21
edm::InputTag
Definition: InputTag.h:15
zPdfUnc::h_thetaZMC_
TH1F * h_thetaZMC_
Definition: zPdfUnc.cc:20
weight
Definition: weight.py:1
zPdfUnc::h_mZMC_
TH1F * h_mZMC_
Definition: zPdfUnc.cc:20
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
zPdfUnc::zPdfUnc
zPdfUnc(const edm::ParameterSet &pset)
Definition: zPdfUnc.cc:49
zPdfUnc::accPtMin_
double accPtMin_
Definition: zPdfUnc.cc:18
Generated for CMSSW Reference Manual by doxygen 1.8.16