File indexing completed on 2024-04-06 12:13:30
0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009
0010 #include "DataFormats/Common/interface/View.h"
0011 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0012
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0016
0017 #include "TH1D.h"
0018
0019 class ExhumeAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0020 public:
0021 explicit ExhumeAnalyzer(const edm::ParameterSet&);
0022 ~ExhumeAnalyzer() override = default;
0023
0024 void analyze(const edm::Event& event, const edm::EventSetup& eventSetup) override;
0025
0026
0027 void beginJob() override;
0028 void endJob() override;
0029
0030 private:
0031 const edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken_;
0032 const double Ebeam_;
0033
0034
0035 TH1D* hist_eta_;
0036 TH1D* hist_phi1_;
0037 TH1D* hist_t1_;
0038 TH1D* hist_xigen1_;
0039 TH1D* hist_phi2_;
0040 TH1D* hist_t2_;
0041 TH1D* hist_xigen2_;
0042 TH1D* hist_sHat_;
0043 TH1D* hist_MX_;
0044 };
0045
0046 ExhumeAnalyzer::ExhumeAnalyzer(const edm::ParameterSet& pset)
0047 : genParticlesToken_(consumes<reco::GenParticleCollection>(pset.getParameter<edm::InputTag>("GenParticleTag"))),
0048 Ebeam_(pset.getParameter<double>("EBeam")) {}
0049
0050 void ExhumeAnalyzer::beginJob() {
0051 edm::Service<TFileService> fs;
0052
0053 hist_eta_ = fs->make<TH1D>("hist_eta", "#eta system", 100, -4.5, 4.5);
0054 hist_phi1_ = fs->make<TH1D>("hist_phi1", "#phi proton 1", 100, -1.1 * M_PI, 1.1 * M_PI);
0055 hist_t1_ = fs->make<TH1D>("hist_t1", "t proton 1", 100, -1.4, 0);
0056 hist_xigen1_ = fs->make<TH1D>("hist_xigen1", "#xi proton 1", 100, 0., 0.1);
0057 hist_phi2_ = fs->make<TH1D>("hist_phi2", "#phi proton 2", 100, -1.1 * M_PI, 1.1 * M_PI);
0058 hist_t2_ = fs->make<TH1D>("hist_t2", "t proton 1", 100, -1.4, 0);
0059 hist_xigen2_ = fs->make<TH1D>("hist_xigen2", "#xi proton 2", 100, 0., 0.1);
0060 hist_sHat_ = fs->make<TH1D>("hist_sHat", "Central inv. mass", 100, 80., 150.);
0061 hist_MX_ = fs->make<TH1D>("hist_MX", "Missing mass", 100, 80., 150.);
0062 }
0063
0064 void ExhumeAnalyzer::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0065
0066 const edm::Handle<reco::GenParticleCollection>& genParticles = event.getHandle(genParticlesToken_);
0067
0068
0069 double pz1max = 0.;
0070 double pz2min = 0.;
0071 reco::GenParticleCollection::const_iterator proton1 = genParticles->end();
0072 reco::GenParticleCollection::const_iterator proton2 = genParticles->end();
0073 for (reco::GenParticleCollection::const_iterator genpart = genParticles->begin(); genpart != genParticles->end();
0074 ++genpart) {
0075 if (genpart->status() != 1)
0076 continue;
0077
0078 double pz = genpart->pz();
0079 if ((genpart->pdgId() == 2212) && (pz > 0.75 * Ebeam_)) {
0080 if (pz > pz1max) {
0081 proton1 = genpart;
0082 pz1max = pz;
0083 }
0084 } else if ((genpart->pdgId() == 2212) && (pz < -0.75 * Ebeam_)) {
0085 if (pz < pz2min) {
0086 proton2 = genpart;
0087 pz2min = pz;
0088 }
0089 }
0090 }
0091
0092 if ((proton1 != genParticles->end()) && (proton2 != genParticles->end())) {
0093 std::cout << "Proton 1: " << proton1->pt() << " " << proton1->eta() << " " << proton1->phi() << std::endl;
0094 std::cout << "Proton 2: " << proton2->pt() << " " << proton2->eta() << " " << proton2->phi() << std::endl;
0095
0096 math::XYZTLorentzVector proton1in(0., 0., Ebeam_, Ebeam_);
0097 math::XYZTLorentzVector proton1diff(proton1->px(), proton1->py(), proton1->pz(), proton1->energy());
0098 double t1 = (proton1diff - proton1in).M2();
0099 double xigen1 = 1 - proton1diff.pz() / Ebeam_;
0100 math::XYZTLorentzVector proton2in(0., 0., -Ebeam_, Ebeam_);
0101 math::XYZTLorentzVector proton2diff(proton2->px(), proton2->py(), proton2->pz(), proton2->energy());
0102 double t2 = (proton2diff - proton2in).M2();
0103 double xigen2 = 1 + proton2diff.pz() / Ebeam_;
0104
0105 double eta = 0.5 * log(xigen1 / xigen2);
0106 double pt1 = sqrt(-t1);
0107 double pt2 = sqrt(-t2);
0108 double phi1 = proton1diff.phi();
0109 double phi2 = proton2diff.phi();
0110 double s = 4 * Ebeam_ * Ebeam_;
0111 double sHat = t1 + t2 - 2 * pt1 * pt2 * cos(phi1 - phi2) + s * xigen1 * xigen2;
0112 double MX = sqrt(s * xigen1 * xigen2);
0113
0114
0115 hist_eta_->Fill(eta);
0116 hist_phi1_->Fill(phi1);
0117 hist_t1_->Fill(t1);
0118 hist_xigen1_->Fill(xigen1);
0119 hist_phi2_->Fill(phi2);
0120 hist_t2_->Fill(t2);
0121 hist_xigen2_->Fill(xigen2);
0122 hist_sHat_->Fill(sqrt(sHat));
0123 hist_MX_->Fill(MX);
0124 }
0125 }
0126
0127 void ExhumeAnalyzer::endJob() {}
0128
0129 #include "FWCore/Framework/interface/MakerMacros.h"
0130 DEFINE_FWK_MODULE(ExhumeAnalyzer);