EventGenerator/interface/HiggsValidation.h

0001 #ifndef HiggsValidation_H
0002 #define HiggsValidation_H
0003
0004 /*class HiggsValidation
0005  *
0006  *  Class to fill Event Generator dqm monitor elements; works on HepMCProduct
0007  *
0008  *
0009  */
0010 #include <iostream>
0011
0012 // framework & common header files
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/Run.h"
0017
0018 #include "DataFormats/Common/interface/Handle.h"
0019 #include "FWCore/Framework/interface/ESHandle.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022
0023 //DQM services
0024 #include "DQMServices/Core/interface/DQMStore.h"
0025 #include "FWCore/ServiceRegistry/interface/Service.h"
0026 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0027
0028 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0029 //#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0030
0031 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0032 #include "TLorentzVector.h"
0033
0034 #include "Validation/EventGenerator/interface/WeightManager.h"
0035
0036 class HiggsValidation : public DQMEDAnalyzer {
0037 public:
0038   explicit HiggsValidation(const edm::ParameterSet &);
0039   ~HiggsValidation() override;
0040
0041   void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override;
0042   void dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) override;
0043   void analyze(edm::Event const &, edm::EventSetup const &) override;
0044
0045 private:
0046   WeightManager wmanager_;
0047
0048   class MonitoredDecays {
0049   public:
0050     MonitoredDecays(const edm::ParameterSet &iConfig) {
0051       fillMap();
0052       std::vector<std::string> input = iConfig.getParameter<std::vector<std::string> >("monitorDecays");
0053       for (std::vector<std::string>::const_iterator i = input.begin(); i != input.end(); ++i) {
0054         fill(*i);
0055       }
0056     }
0057
0058     ~MonitoredDecays(){};
0059
0060     size_t position(int pid1, int pid2) {
0061       if (abs(pid1) == 14 || abs(pid1) == 16)
0062         pid1 = 12;
0063       if (abs(pid2) == 14 || abs(pid2) == 16)
0064         pid2 = 12;
0065       for (size_t i = 0; i < channels.size(); ++i) {
0066         if ((channels[i].first == abs(pid1) && channels[i].second == abs(pid2)) ||
0067             (channels[i].first == abs(pid2) && channels[i].second == abs(pid1)))
0068           return i + 1;
0069       }
0070       return undetermined();  //channels.size()+1;
0071     }
0072
0073     size_t size() { return channels.size() + 2; }
0074     size_t undetermined() { return 0; }
0075     size_t stable() { return size(); }
0076     std::string channel(size_t i) {
0077       if (i == 0)
0078         return "?";
0079       if (i == channels.size() + 1)
0080         return "Stable";
0081       return convert(channels[i - 1].first) + convert(channels[i - 1].second);
0082     }
0083
0084     int convert(std::string s) {
0085       if (namePidMap.count(s)) {
0086         return namePidMap[s];
0087       }
0088       return 0;
0089     }
0090
0091     std::string convert(int pid) {
0092       pid = abs(pid);
0093       if (pid == 14 || pid == 16)
0094         pid = 12;
0095       for (std::map<std::string, int>::const_iterator i = namePidMap.begin(); i != namePidMap.end(); ++i) {
0096         if (i->second == pid)
0097           return i->first;
0098       }
0099       return "not found";
0100     }
0101
0102     unsigned int NDecayParticles() { return nparticles_; }
0103
0104     int isDecayParticle(int pid) {
0105       int idx = 0;
0106       for (std::map<std::string, int>::const_iterator i = namePidMap.begin(); i != namePidMap.end(); ++i) {
0107         if (i->second == pid)
0108           return idx;
0109         idx++;
0110       }
0111       return -1;
0112     }
0113
0114     std::string ConvertIndex(int index) {
0115       int idx = 0;
0116       for (std::map<std::string, int>::const_iterator i = namePidMap.begin(); i != namePidMap.end(); ++i) {
0117         if (idx == index)
0118           return i->first;
0119         idx++;
0120       }
0121       return "unknown";
0122     }
0123
0124   private:
0125     void fill(std::string s) {
0126       size_t pos = s.find('+');
0127       std::string particle1 = s.substr(0, pos);
0128       std::string particle2 = s.substr(pos + 1, s.length() - pos);
0129       std::pair<int, int> decay;
0130       decay.first = convert(particle1);
0131       decay.second = convert(particle2);
0132       channels.push_back(decay);
0133     }
0134
0135     void fillMap() {
0136       namePidMap["d"] = 1;
0137       namePidMap["u"] = 2;
0138       namePidMap["s"] = 3;
0139       namePidMap["c"] = 4;
0140       namePidMap["b"] = 5;
0141       namePidMap["t"] = 6;
0142       namePidMap["e"] = 11;
0143       namePidMap["nu"] = 12;
0144       namePidMap["mu"] = 13;
0145       namePidMap["tau"] = 15;
0146       namePidMap["gamma"] = 22;
0147       namePidMap["Z"] = 23;
0148       namePidMap["W"] = 24;
0149       nparticles_ = 0;
0150       for (std::map<std::string, int>::const_iterator i = namePidMap.begin(); i != namePidMap.end(); ++i) {
0151         nparticles_++;
0152       }
0153     }
0154
0155     std::map<std::string, int> namePidMap;
0156
0157     std::vector<std::pair<int, int> > channels;
0158     unsigned int nparticles_;
0159   };
0160
0161   int findHiggsDecayChannel(const HepMC::GenParticle *, std::vector<HepMC::GenParticle *> &decayprod);
0162   std::string convert(int);
0163
0164   edm::InputTag hepmcCollection_;
0165
0166   int particle_id;
0167   std::string particle_name;
0168
0169   MonitoredDecays *monitoredDecays;
0170
0171   /// PDT table
0172   edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable;
0173   edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> fPDGTableToken;
0174
0175   MonitorElement *nEvt;
0176   MonitorElement *HiggsDecayChannels;
0177
0178   MonitorElement *Higgs_pt;
0179   MonitorElement *Higgs_eta;
0180   MonitorElement *Higgs_mass;
0181
0182   std::vector<MonitorElement *> HiggsDecayProd_pt;
0183   std::vector<MonitorElement *> HiggsDecayProd_eta;
0184
0185   edm::EDGetTokenT<edm::HepMCProduct> hepmcCollectionToken_;
0186 };
0187
0188 #endif