BasicHepMCValidation.h

Go to the documentation of this file.

#ifndef BASICHEPMCVALIDATION_H
#define BASICHEPMCVALIDATION_H

/*class BasicHepMCValidation
 *  
 *  Class to fill Event Generator dqm monitor elements; works on HepMCProduct
 *
 *
 */

// framework & common header files
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Run.h"

#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"

//DQM services
#include "DQMServices/Core/interface/DQMStore.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
#include "Validation/EventGenerator/interface/DQMHelper.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"

#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"

#include "Validation/EventGenerator/interface/WeightManager.h"
#include "TVector3.h"

class BasicHepMCValidation : public DQMEDAnalyzer {
public:
  explicit BasicHepMCValidation(const edm::ParameterSet &);
  ~BasicHepMCValidation() override;

  void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override;
  void dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) override;
  void analyze(edm::Event const &, edm::EventSetup const &) override;

private:
  WeightManager wmanager_;
  edm::InputTag hepmcCollection_;

  edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable;
  edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> fPDGTableToken;

  class ParticleMonitor {
  public:
    ParticleMonitor(std::string name_, int pdgid_, DQMStore::IBooker &i, bool nlog_ = false)
        : name(name_), pdgid(pdgid_), count(0), nlog(nlog_) {
      DQMHelper dqm(&i);
      // Number of analyzed events
      if (!nlog) {
        numberPerEvent = dqm.book1dHisto(
            name + "Number", "Number of  " + name + "'s per event", 20, 0, 20, "No. of " + name, "Number of Events");
      } else {
        numberPerEvent = dqm.book1dHisto(name + "Number",
                                         "Number of  " + name + "'s per event",
                                         20,
                                         0,
                                         20,
                                         "log_{10}(No. of " + name + ")",
                                         "Number of Events");
      }
      p_init = dqm.book1dHisto(name + "Momentum",
                               "log_{10}(P) of the " + name + "s",
                               60,
                               -2,
                               4,
                               "log_{10}(P) (log_{10}(GeV))",
                               "Number of " + name);

      eta_init = dqm.book1dHisto(name + "Eta", "#eta of the " + name + "s", 100, -5., 5., "#eta", "Number of " + name);

      lifetime_init = dqm.book1dHisto(name + "LifeTime",
                                      "#phi of the " + name + "s",
                                      100,
                                      -15,
                                      -5,
                                      "Log_{10}(life-time^{final}) (log_{10}(s))",
                                      "Number of " + name);

      p_final = dqm.book1dHisto(name + "MomentumFinal",
                                "log_{10}(P^{final}) of " + name + "s at end of decay chain",
                                60,
                                -2,
                                4,
                                "log_{10}(P^{final}) (log_{10}(GeV))",
                                "Number of " + name);

      lifetime_final = dqm.book1dHisto(name + "LifeTimeFinal",
                                       "Log_{10}(life-time^{final}) of " + name + "s at end of decay chain",
                                       100,
                                       -15,
                                       -5,
                                       "Log_{10}(life-time^{final}) (log_{10}(s))",
                                       "Number of " + name);
    }

    ~ParticleMonitor(){};

    bool Fill(const HepMC::GenParticle *p, double weight) {
      if (p->pdg_id() == pdgid) {
        if (isFirst(p)) {
          p_init->Fill(log10(p->momentum().rho()), weight);
          eta_init->Fill(p->momentum().eta(), weight);
          const HepMC::GenParticle *pf = GetFinal(p);  // inlcude mixing
          p_final->Fill(log10(pf->momentum().rho()), weight);
          // compute lifetime...
          if (p->production_vertex() && p->end_vertex()) {
            TVector3 PV(p->production_vertex()->point3d().x(),
                        p->production_vertex()->point3d().y(),
                        p->production_vertex()->point3d().z());
            TVector3 SV(p->end_vertex()->point3d().x(), p->end_vertex()->point3d().y(), p->end_vertex()->point3d().z());
            TVector3 DL = SV - PV;
            double c(2.99792458E8), Ltau(DL.Mag() / 100) /*cm->m*/, beta(p->momentum().rho() / p->momentum().m());
            double lt = Ltau / (c * beta);
            if (lt > 1E-16)
              lifetime_init->Fill(log10(lt), weight);
            if (pf->end_vertex()) {
              TVector3 SVf(
                  pf->end_vertex()->point3d().x(), pf->end_vertex()->point3d().y(), pf->end_vertex()->point3d().z());
              DL = SVf - PV;
              Ltau = DL.Mag() / 100;
              lt = Ltau / (c * beta);
              if (lt > 1E-16)
                lifetime_final->Fill(log10(lt), weight);
            }
          }
          count++;
        }
        return true;
      }
      return false;
    }

    void FillCount(double weight) {
      if (nlog)
        numberPerEvent->Fill(log10(count), weight);
      else
        numberPerEvent->Fill(count, weight);
      count = 0;
    }

    int PDGID() { return pdgid; }

  private:
    bool isFirst(const HepMC::GenParticle *p) {
      if (p->production_vertex()) {
        for (HepMC::GenVertex::particles_in_const_iterator m = p->production_vertex()->particles_in_const_begin();
             m != p->production_vertex()->particles_in_const_end();
             m++) {
          if (abs((*m)->pdg_id()) == abs(p->pdg_id()))
            return false;
        }
      }
      return true;
    }

    // includes mixing (assuming mixing is not occurring more than 5 times back and forth)
    const HepMC::GenParticle *GetFinal(const HepMC::GenParticle *p) {
      const HepMC::GenParticle *aPart = p;
      for (unsigned int iMix = 0; iMix < 10; iMix++) {
        bool foundSimilar = false;
        if (aPart->end_vertex()) {
          if (aPart->end_vertex()->particles_out_size() != 0) {
            for (HepMC::GenVertex::particles_out_const_iterator d = aPart->end_vertex()->particles_out_const_begin();
                 d != aPart->end_vertex()->particles_out_const_end();
                 d++) {
              if (abs((*d)->pdg_id()) == abs(aPart->pdg_id())) {
                aPart = *d;
                foundSimilar = true;
                break;
              }
            }
          }
          if (!foundSimilar)
            break;
        }
      }
      return aPart;
    }

    std::string name;
    int pdgid;
    unsigned int count;
    bool nlog;
    MonitorElement *p_init, *p_final, *eta_init, *lifetime_init, *lifetime_final, *numberPerEvent;
  };

  MonitorElement *nEvt;
  std::vector<ParticleMonitor> particles;

  MonitorElement *otherPtclNumber;
  MonitorElement *otherPtclMomentum;
  MonitorElement *genPtclNumber;
  MonitorElement *genVrtxNumber;
  MonitorElement *unknownPDTNumber;
  MonitorElement *outVrtxPtclNumber;
  MonitorElement *genPtclStatus;
  //
  MonitorElement *stablePtclNumber;
  MonitorElement *stableChaNumber;
  MonitorElement *stablePtclPhi;
  MonitorElement *stablePtclEta;
  MonitorElement *stablePtclCharge;
  MonitorElement *stablePtclp;
  MonitorElement *stablePtclpT;
  MonitorElement *partonNumber;
  MonitorElement *partonpT;
  MonitorElement *outVrtxStablePtclNumber;
  //
  MonitorElement *vrtxZ;
  MonitorElement *vrtxRadius;
  //
  MonitorElement *Bjorken_x;
  MonitorElement *pdf_u;
  MonitorElement *pdf_ubar;
  MonitorElement *pdf_d;
  MonitorElement *pdf_dbar;
  MonitorElement *pdf_ssbar;
  MonitorElement *pdf_ccbar;
  MonitorElement *pdf_bbbar;
  MonitorElement *pdf_g;
  MonitorElement *scalePDF;
  MonitorElement *parton1Id;
  MonitorElement *parton2Id;

  MonitorElement *status1ShortLived;

  MonitorElement *log10DeltaEcms;
  MonitorElement *DeltaEcms;
  MonitorElement *DeltaPx;
  MonitorElement *DeltaPy;
  MonitorElement *DeltaPz;

  edm::EDGetTokenT<edm::HepMCProduct> hepmcCollectionToken_;
};

#endif