ExhumeInterface/test/ExhumeAnalyzer.cc

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   //virtual void beginJob(const edm::EventSetup& eventSetup);
0027   void beginJob() override;
0028   void endJob() override;
0029
0030 private:
0031   const edm::EDGetTokenT<reco::GenParticleCollection> genParticlesToken_;
0032   const double Ebeam_;
0033
0034   // Histograms
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   // Generator Information
0066   const edm::Handle<reco::GenParticleCollection>& genParticles = event.getHandle(genParticlesToken_);
0067
0068   // Look for protons
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     //Fill histograms
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);