File indexing completed on 2024-04-06 12:13:55
0001 #include <iostream>
0002
0003 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0004 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0005 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0006
0007
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "DataFormats/Common/interface/Handle.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0014 #include "TH1.h"
0015
0016 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0017
0018 class HZZ4muAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0019 public:
0020
0021 explicit HZZ4muAnalyzer(const edm::ParameterSet&);
0022 ~HZZ4muAnalyzer() = default;
0023
0024
0025 void analyze(const edm::Event&, const edm::EventSetup&) override;
0026 void beginJob() override;
0027
0028 private:
0029 const edm::EDGetTokenT<GenEventInfoProduct> tokenGenInfo_;
0030 const edm::EDGetTokenT<edm::HepMCProduct> tokenHepMC_;
0031 TH1D* fHist2muMass;
0032 TH1D* fHist4muMass;
0033 TH1D* fHistZZMass;
0034 };
0035
0036 HZZ4muAnalyzer::HZZ4muAnalyzer(const edm::ParameterSet& pset)
0037 : tokenGenInfo_(consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
0038 tokenHepMC_(consumes<edm::HepMCProduct>(edm::InputTag("VtxSmeared"))),
0039 fHist2muMass(nullptr),
0040 fHist4muMass(nullptr),
0041 fHistZZMass(nullptr) {
0042
0043 usesResource(TFileService::kSharedResource);
0044 }
0045
0046 void HZZ4muAnalyzer::beginJob() {
0047 edm::Service<TFileService> fs;
0048 fHist2muMass = fs->make<TH1D>("Hist2muMass", "2-mu inv. mass", 100, 60., 120.);
0049 fHist4muMass = fs->make<TH1D>("Hist4muMass", "4-mu inv. mass", 100, 170., 210.);
0050 fHistZZMass = fs->make<TH1D>("HistZZMass", "ZZ inv. mass", 100, 170., 210.);
0051 }
0052
0053 void HZZ4muAnalyzer::analyze(const edm::Event& e, const edm::EventSetup&) {
0054
0055 const edm::Handle<GenEventInfoProduct>& GenInfoHandle = e.getHandle(tokenGenInfo_);
0056 double qScale = GenInfoHandle->qScale();
0057 double pthat = (GenInfoHandle->hasBinningValues() ? (GenInfoHandle->binningValues())[0] : 0.0);
0058 std::cout << " qScale = " << qScale << " pthat = " << pthat << std::endl;
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075 const edm::Handle<edm::HepMCProduct>& EvtHandle = e.getHandle(tokenHepMC_);
0076
0077 const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0078
0079
0080
0081
0082
0083 HepMC::GenVertex* HiggsDecVtx = 0;
0084
0085
0086
0087
0088
0089
0090
0091 for (HepMC::GenEvent::vertex_const_iterator vit = Evt->vertices_begin(); vit != Evt->vertices_end(); vit++) {
0092 for (HepMC::GenVertex::particles_out_const_iterator pout = (*vit)->particles_out_const_begin();
0093 pout != (*vit)->particles_out_const_end();
0094 pout++) {
0095 if ((*pout)->pdg_id() == 25) {
0096 if ((*pout)->end_vertex() != 0) {
0097 HiggsDecVtx = (*pout)->end_vertex();
0098 break;
0099 }
0100 }
0101 }
0102 if (HiggsDecVtx != 0) {
0103 break;
0104 }
0105 }
0106
0107 if (HiggsDecVtx == 0) {
0108 std::cout << " There is NO Higgs in this event ! " << std::endl;
0109 return;
0110 }
0111
0112 if (e.id().event() == 1) {
0113 std::cout << " " << std::endl;
0114 std::cout << " We do some example printouts in the event 1 " << std::endl;
0115 std::cout << " Higgs decay found at the vertex " << HiggsDecVtx->barcode() << " (barcode)" << std::endl;
0116
0117 std::vector<HepMC::GenParticle*> HiggsChildren;
0118
0119 for (HepMC::GenVertex::particles_out_const_iterator H0in = HiggsDecVtx->particles_out_const_begin();
0120 H0in != HiggsDecVtx->particles_out_const_end();
0121 H0in++) {
0122 HiggsChildren.push_back(*H0in);
0123 }
0124 std::cout << " Number of Higgs (immediate) children = " << HiggsChildren.size() << std::endl;
0125 for (unsigned int ic = 0; ic < HiggsChildren.size(); ic++) {
0126 HiggsChildren[ic]->print();
0127 }
0128 }
0129
0130
0131
0132 std::vector<HepMC::GenParticle*> StableHiggsDesc;
0133
0134 if (e.id().event() == 1)
0135 std::cout << " Now let us list all descendents of the Higgs" << std::endl;
0136 for (HepMC::GenVertex::particle_iterator des = HiggsDecVtx->particles_begin(HepMC::descendants);
0137 des != HiggsDecVtx->particles_end(HepMC::descendants);
0138 des++) {
0139 if (e.id().event() == 1)
0140 (*des)->print();
0141 if ((*des)->status() == 1)
0142 StableHiggsDesc.push_back(*des);
0143 }
0144
0145 HepMC::FourVector Mom2part;
0146 double XMass2part = 0.;
0147 double XMass4part = 0.;
0148 double XMass2pairs = 0.;
0149 std::vector<HepMC::FourVector> Mom2partCont;
0150
0151
0152
0153
0154 for (unsigned int i = 0; i < StableHiggsDesc.size(); i++) {
0155
0156
0157 if (std::abs(StableHiggsDesc[i]->pdg_id()) != 13)
0158 continue;
0159
0160 for (unsigned int j = i + 1; j < StableHiggsDesc.size(); j++) {
0161
0162
0163 if (std::abs(StableHiggsDesc[j]->pdg_id()) != 13)
0164 continue;
0165
0166
0167
0168 if ((StableHiggsDesc[i]->pdg_id() * StableHiggsDesc[j]->pdg_id()) > 0)
0169 continue;
0170
0171
0172
0173 Mom2part = HepMC::FourVector((StableHiggsDesc[i]->momentum().px() + StableHiggsDesc[j]->momentum().px()),
0174 (StableHiggsDesc[i]->momentum().py() + StableHiggsDesc[j]->momentum().py()),
0175 (StableHiggsDesc[i]->momentum().pz() + StableHiggsDesc[j]->momentum().pz()),
0176 (StableHiggsDesc[i]->momentum().e() + StableHiggsDesc[j]->momentum().e()));
0177
0178 XMass2part = Mom2part.m();
0179 fHist2muMass->Fill(XMass2part);
0180
0181
0182
0183
0184
0185
0186 if (XMass2part > 80. && XMass2part < 100.) {
0187 Mom2partCont.push_back(Mom2part);
0188 }
0189 }
0190 }
0191
0192
0193
0194
0195 double px4, py4, pz4, e4;
0196 px4 = py4 = pz4 = e4 = 0.;
0197 if (StableHiggsDesc.size() == 4) {
0198 for (unsigned int i = 0; i < StableHiggsDesc.size(); i++) {
0199 px4 += StableHiggsDesc[i]->momentum().px();
0200 py4 += StableHiggsDesc[i]->momentum().py();
0201 pz4 += StableHiggsDesc[i]->momentum().pz();
0202 e4 += StableHiggsDesc[i]->momentum().e();
0203 }
0204 XMass4part = HepMC::FourVector(px4, py4, pz4, e4).m();
0205 fHist4muMass->Fill(XMass4part);
0206 }
0207
0208
0209
0210
0211
0212 for (unsigned int i = 0; i < Mom2partCont.size(); i++) {
0213 for (unsigned int j = i + 1; j < Mom2partCont.size(); j++) {
0214
0215 XMass2pairs = HepMC::FourVector((Mom2partCont[i].px() + Mom2partCont[j].px()),
0216 (Mom2partCont[i].py() + Mom2partCont[j].py()),
0217 (Mom2partCont[i].pz() + Mom2partCont[j].pz()),
0218 (Mom2partCont[i].e() + Mom2partCont[j].e()))
0219 .m();
0220 fHistZZMass->Fill(XMass2pairs);
0221
0222 }
0223 }
0224
0225 return;
0226 }
0227
0228 typedef HZZ4muAnalyzer HZZ4muExampleAnalyzer;
0229 DEFINE_FWK_MODULE(HZZ4muExampleAnalyzer);