BasicHepMCValidation.cc

Go to the documentation of this file.

/*class BasicHepMCValidation
 *  
 *  Class to fill dqm monitor elements from existing EDM file
 *
 */

#include "Validation/EventGenerator/interface/BasicHepMCValidation.h"

#include "CLHEP/Units/defs.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include "Validation/EventGenerator/interface/DQMHelper.h"
using namespace edm;

BasicHepMCValidation::BasicHepMCValidation(const edm::ParameterSet &iPSet)
    : wmanager_(iPSet, consumesCollector()), hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection")) {
  hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_);
}

BasicHepMCValidation::~BasicHepMCValidation() {}

void BasicHepMCValidation::dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) { c.getData(fPDGTable); }

// Set upper bound & lower bound for PDF & Scale related histograms
double logPdfMax = 3.0;
double logPdfMin = -3.0;
int logPdfNbin = 150;
double logPdfBinsize = (logPdfMax - logPdfMin) / (double)logPdfNbin;
double logQScaleMax = 4.0;
double logQScaleMin = -1.0;
int logQScaleNbin = 500;
double logQScaleBinsize = (logQScaleMax - logQScaleMin) / (double)logQScaleNbin;

void BasicHepMCValidation::bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) {
  DQMHelper dqm(&i);
  i.setCurrentFolder("Generator/Particles");

  // Number of analyzed events
  nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1.);

  // quarks
  particles.push_back(ParticleMonitor("u", 1, i));
  particles.push_back(ParticleMonitor("ubar", -1, i));
  particles.push_back(ParticleMonitor("d", 2, i));
  particles.push_back(ParticleMonitor("dbar", -2, i));
  particles.push_back(ParticleMonitor("s", 3, i));
  particles.push_back(ParticleMonitor("sbar", -3, i));
  particles.push_back(ParticleMonitor("c", 4, i));
  particles.push_back(ParticleMonitor("cbar", -4, i));
  particles.push_back(ParticleMonitor("b", 5, i));
  particles.push_back(ParticleMonitor("bbar", -5, i));
  particles.push_back(ParticleMonitor("t", 6, i));
  particles.push_back(ParticleMonitor("tbar", -6, i));

  //leptons
  particles.push_back(ParticleMonitor("eminus", 11, i));
  particles.push_back(ParticleMonitor("eplus", -11, i));
  particles.push_back(ParticleMonitor("nue", 12, i));
  particles.push_back(ParticleMonitor("nuebar", -12, i));
  particles.push_back(ParticleMonitor("muminus", 13, i));
  particles.push_back(ParticleMonitor("muplus", -13, i));
  particles.push_back(ParticleMonitor("numu", 14, i));
  particles.push_back(ParticleMonitor("numubar", -14, i));
  particles.push_back(ParticleMonitor("tauminus", 15, i));
  particles.push_back(ParticleMonitor("tauplus", -15, i));
  particles.push_back(ParticleMonitor("nutau", 16, i));
  particles.push_back(ParticleMonitor("nutaubar", -16, i));

  //bosons
  particles.push_back(ParticleMonitor("Wplus", 24, i));
  particles.push_back(ParticleMonitor("Wminus", -24, i));
  particles.push_back(ParticleMonitor("Z", 23, i));
  particles.push_back(ParticleMonitor("gamma", 22, i));
  particles.push_back(ParticleMonitor("gluon", 21, i));

  //mesons
  particles.push_back(ParticleMonitor("piplus", 211, i, true));    //log
  particles.push_back(ParticleMonitor("piminus", -211, i, true));  //log
  particles.push_back(ParticleMonitor("pizero", 111, i, true));    //log
  particles.push_back(ParticleMonitor("Kplus", 321, i));
  particles.push_back(ParticleMonitor("Kminus", -321, i));
  particles.push_back(ParticleMonitor("Klzero", 130, i));
  particles.push_back(ParticleMonitor("Kszero", 310, i));

  //baryons
  particles.push_back(ParticleMonitor("p", 2212, i, true));      //log
  particles.push_back(ParticleMonitor("pbar", -2212, i, true));  //log
  particles.push_back(ParticleMonitor("n", 2112, i, true));      //log
  particles.push_back(ParticleMonitor("nbar", -2112, i, true));  //log
  particles.push_back(ParticleMonitor("lambda0", 3122, i));
  particles.push_back(ParticleMonitor("lambda0bar", -3122, i));

  //D mesons
  particles.push_back(ParticleMonitor("Dplus", 411, i));
  particles.push_back(ParticleMonitor("Dminus", -411, i));
  particles.push_back(ParticleMonitor("Dzero", 421, i));
  particles.push_back(ParticleMonitor("Dzerobar", -421, i));

  //B mesons
  particles.push_back(ParticleMonitor("Bplus", 521, i));
  particles.push_back(ParticleMonitor("Bminus", -521, i));
  particles.push_back(ParticleMonitor("Bzero", 511, i));
  particles.push_back(ParticleMonitor("Bzerobar", -511, i));
  particles.push_back(ParticleMonitor("Bszero", 531, i));
  particles.push_back(ParticleMonitor("Bszerobar", -531, i));

  //
  otherPtclNumber = dqm.book1dHisto(
      "otherPtclNumber", "Log10(No. other ptcls)", 60, -1, 5, "log_{10}(No. other ptcls)", "Number of Events");  //Log
  otherPtclMomentum = dqm.book1dHisto("otherPtclMomentum",
                                      "Log10(p) other ptcls",
                                      60,
                                      -2,
                                      4,
                                      "log10(P^{other ptcls}) (log_{10}(GeV))",
                                      "Number of Events");

  genPtclNumber = dqm.book1dHisto(
      "genPtclNumber", "Log10(No. all particles)", 60, -1, 5, "log10(No. all particles)", "Number of Events");  //Log
  genVrtxNumber = dqm.book1dHisto(
      "genVrtxNumber", "Log10(No. all vertexs)", 60, -1, 5, "log10(No. all vertexs)", "Number of Events");  //Log
  //
  stablePtclNumber = dqm.book1dHisto("stablePtclNumber",
                                     "Log10(No. stable particles)",
                                     50,
                                     0,
                                     5,
                                     "log10(No. stable particles)",
                                     "Number of Events");  //Log
  stablePtclPhi = dqm.book1dHisto(
      "stablePtclPhi", "stable Ptcl Phi", 360, -180, 180, "#phi^{stable Ptcl} (rad)", "Number of Events");
  stablePtclEta = dqm.book1dHisto(
      "stablePtclEta", "stable Ptcl Eta (pseudo rapidity)", 220, -11, 11, "#eta^{stable Ptcl}", "Number of Events");
  stablePtclCharge =
      dqm.book1dHisto("stablePtclCharge", "stablePtclCharge", 5, -2, 2, "Charge^{stable ptcls}", "Number of Events");
  stableChaNumber = dqm.book1dHisto("stableChaNumber",
                                    "Log10(No. stable charged particles)",
                                    50,
                                    0,
                                    5,
                                    "log_{10}(No. stable charged particles)",
                                    "Number of Events");  //Log
  stablePtclp = dqm.book1dHisto("stablePtclp",
                                "Log10(p) stable ptcl p",
                                80,
                                -4,
                                4,
                                "log_{10}(P^{stable ptcl}) (log_{10}(GeV))",
                                "Number of Events");  //Log
  stablePtclpT = dqm.book1dHisto("stablePtclpT",
                                 "Log10(pT) stable ptcl pT",
                                 80,
                                 -4,
                                 4,
                                 "log_{10}(P_{t}^{stable ptcl}) (log_{10}(GeV))",
                                 "Number of Events");  //Log
  partonNumber =
      dqm.book1dHisto("partonNumber", "number of partons", 100, 0, 100, "number of partons", "Number of Events");
  partonpT =
      dqm.book1dHisto("partonpT", "Log10(pT) parton pT", 80, -4, 4, "Log10(P_{t}^{parton})", "Number of Events");  //Log
  outVrtxStablePtclNumber = dqm.book1dHisto("outVrtxStablePtclNumber",
                                            "No. outgoing stable ptcls from vrtx",
                                            10,
                                            0,
                                            10,
                                            "No. outgoing stable ptcls from vrtx",
                                            "Number of Events");
  //
  outVrtxPtclNumber = dqm.book1dHisto("outVrtxPtclNumber",
                                      "No. outgoing ptcls from vrtx",
                                      30,
                                      0,
                                      30,
                                      "No. outgoing ptcls from vrtx",
                                      "Number of Events");
  vrtxZ = dqm.book1dHisto("VrtxZ", "VrtxZ", 50, -250, 250, "Z_{Vtx}", "Number of Events");
  vrtxRadius = dqm.book1dHisto("vrtxRadius", "vrtxRadius", 50, 0, 50, "R_{vtx}", "Number of Events");
  //
  unknownPDTNumber = dqm.book1dHisto("unknownPDTNumber",
                                     "Log10(No. unknown ptcls PDT)",
                                     60,
                                     -1,
                                     5,
                                     "log_{10}(No. unknown ptcls PDT)",
                                     "Number of Events");  //Log
  genPtclStatus = dqm.book1dHisto("genPtclStatus", "Status of genParticle", 200, 0, 200., "", "Number of Events");
  //
  Bjorken_x = dqm.book1dHisto("Bjorken_x", "Bjorken_x", 1000, 0.0, 1.0, "Bjorken_{x}", "Number of Events");
  pdf_u = dqm.book1dHisto(
      "pdf_u", "Log10(PDF(u,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events");  //Log
  pdf_ubar = dqm.book1dHisto("pdf_ubar",
                             "Log10(PDF(ubar,x,Q))",
                             logPdfNbin,
                             logPdfMin,
                             logPdfMax,
                             "log_{10}(x*f(x))",
                             "Number of Events");  //Log
  pdf_d = dqm.book1dHisto(
      "pdf_d", "Log10(PDF(d,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events");  //Log
  pdf_dbar = dqm.book1dHisto("pdf_dbar",
                             "Log10(PDF(dbar,x,Q))",
                             logPdfNbin,
                             logPdfMin,
                             logPdfMax,
                             "log_{10}(x*f(x))",
                             "Number of Events");  //Log
  pdf_ssbar = dqm.book1dHisto("pdf_ssbar",
                              "Log10(PDF(ssbar,x,Q))",
                              logPdfNbin,
                              logPdfMin,
                              logPdfMax,
                              "log_{10}(x*f(x))",
                              "Number of Events");  //Log
  pdf_ccbar = dqm.book1dHisto("pdf_ccbar",
                              "Log10(PDF(ccbar,x,Q))",
                              logPdfNbin,
                              logPdfMin,
                              logPdfMax,
                              "log_{10}(x*f(x))",
                              "Number of Events");  //Log
  pdf_bbbar = dqm.book1dHisto("pdf_bbbar",
                              "Log10(PDF(bbbar,x,Q))",
                              logPdfNbin,
                              logPdfMin,
                              logPdfMax,
                              "log_{10}(x*f(x))",
                              "Number of Events");  //Log
  pdf_g = dqm.book1dHisto(
      "pdf_g", "Log10(PDF(g,x,Q))", logPdfNbin, logPdfMin, logPdfMax, "log_{10}(x*f(x))", "Number of Events");  //Log
  scalePDF = dqm.book1dHisto("scalePDF",
                             "Log10(Q-scale(GeV))",
                             logQScaleNbin,
                             logQScaleMin,
                             logQScaleMax,
                             "log_{10}(Q-scale(GeV))",
                             "Number of Events");
  parton1Id = dqm.book1dHisto("parton1Id", "ID of parton 1", 45, -14.5, 30.5, "ID", "Number of Events");
  parton2Id = dqm.book1dHisto("parton2Id", "ID of parton 2", 45, -14.5, 30.5, "ID", "Number of Events");
  //
  status1ShortLived = dqm.book1dHisto("status1ShortLived", "Status 1 short lived", 11, 0, 11, "", "Number of Events");
  status1ShortLived->setBinLabel(1, "d/dbar");
  status1ShortLived->setBinLabel(2, "u/ubar");
  status1ShortLived->setBinLabel(3, "s/sbar");
  status1ShortLived->setBinLabel(4, "c/cbar");
  status1ShortLived->setBinLabel(5, "b/bbar");
  status1ShortLived->setBinLabel(6, "t/tbar");
  status1ShortLived->setBinLabel(7, "g");
  status1ShortLived->setBinLabel(8, "tau-/tau+");
  status1ShortLived->setBinLabel(9, "Z0");
  status1ShortLived->setBinLabel(10, "W-/W+");
  status1ShortLived->setBinLabel(11, "PDG = 7,8,17,25-99");

  log10DeltaEcms = dqm.book1dHisto("DeltaEcms1log10",
                                   "log_{10} of deviation from nominal Ecms",
                                   200,
                                   -1.,
                                   5.,
                                   "log_{10}(#DeltaE) (log_{10}(GeV))",
                                   "Number of Events");
  DeltaEcms =
      dqm.book1dHisto("DeltaEcms1", "deviation from nominal Ecms", 200, -1., 1., "#DeltaE (GeV)", "Number of Events");
  DeltaPx =
      dqm.book1dHisto("DeltaPx1", "deviation from nominal Px", 200, -1., 1., "#DeltaP_{x} (GeV)", "Number of Events");
  DeltaPy =
      dqm.book1dHisto("DeltaPy1", "deviation from nominal Py", 200, -1., 1., "#DeltaP_{y} (GeV)", "Number of Events");
  DeltaPz =
      dqm.book1dHisto("DeltaPz1", "deviation from nominal Pz", 200, -1., 1., "#DeltaP_{z} (GeV)", "Number of Events");
  return;
}

void BasicHepMCValidation::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
  int partonNum = 0;
  //
  int outVrtxStablePtclNum = 0;
  int stablePtclNum = 0;
  int otherPtclNum = 0;
  int unknownPDTNum = 0;
  int stableChaNum = 0;
  //
  double bjorken = 0.;
  double logPdf1 = 0.;
  double logPdf2 = 0.;
  double logQScale = 0.;
  //
  double etotal = 0.;
  double pxtotal = 0.;
  double pytotal = 0.;
  double pztotal = 0.;

  edm::Handle<HepMCProduct> evt;
  iEvent.getByToken(hepmcCollectionToken_, evt);

  //Get EVENT
  HepMC::GenEvent const *myGenEvent = evt->GetEvent();

  double weight = wmanager_.weight(iEvent);

  nEvt->Fill(0.5, weight);

  genPtclNumber->Fill(log10(myGenEvent->particles_size()), weight);
  genVrtxNumber->Fill(log10(myGenEvent->vertices_size()), weight);

  HepMC::PdfInfo const *pdf = myGenEvent->pdf_info();
  if (pdf) {
    bjorken = ((pdf->x1()) / ((pdf->x1()) + (pdf->x2())));
    logQScale = log10(pdf->scalePDF());
    if (logQScale > logQScaleMax)
      logQScale = logQScaleMax - 0.5 * logQScaleBinsize;  // visualize overflow & underflow in the histograms
    if (logQScale < logQScaleMin)
      logQScale = logQScaleMin + 0.5 * logQScaleBinsize;
    logPdf1 = log10(pdf->pdf1());
    if (logPdf1 > logPdfMax)
      logPdf1 = logPdfMax - 0.5 * logPdfBinsize;
    if (logPdf1 < logPdfMin)
      logPdf1 = logPdfMin + 0.5 * logPdfBinsize;
    logPdf2 = log10(pdf->pdf2());
    if (logPdf2 > logPdfMax)
      logPdf2 = logPdfMax - 0.5 * logPdfBinsize;
    if (logPdf2 < logPdfMin)
      logPdf2 = logPdfMin + 0.5 * logPdfBinsize;
    Bjorken_x->Fill(bjorken, weight);
    scalePDF->Fill(logQScale, weight);
    parton1Id->Fill((double)pdf->id1(), weight);
    parton2Id->Fill((double)pdf->id2(), weight);
    if (pdf->id1() == 2)
      pdf_u->Fill(logPdf1, weight);
    if (pdf->id2() == 2)
      pdf_u->Fill(logPdf2, weight);
    if (pdf->id1() == -2)
      pdf_ubar->Fill(logPdf1, weight);
    if (pdf->id2() == -2)
      pdf_ubar->Fill(logPdf2, weight);
    if (pdf->id1() == 1)
      pdf_d->Fill(logPdf1, weight);
    if (pdf->id2() == 1)
      pdf_d->Fill(logPdf2, weight);
    if (pdf->id1() == -1)
      pdf_dbar->Fill(logPdf1, weight);
    if (pdf->id2() == -1)
      pdf_dbar->Fill(logPdf2, weight);
    if (std::abs(pdf->id1()) == 3)
      pdf_ssbar->Fill(logPdf1, weight);
    if (std::abs(pdf->id2()) == 3)
      pdf_ssbar->Fill(logPdf2, weight);
    if (std::abs(pdf->id1()) == 4)
      pdf_ccbar->Fill(logPdf1, weight);
    if (std::abs(pdf->id2()) == 4)
      pdf_ccbar->Fill(logPdf2, weight);
    if (std::abs(pdf->id1()) == 5)
      pdf_bbbar->Fill(logPdf1, weight);
    if (std::abs(pdf->id2()) == 5)
      pdf_bbbar->Fill(logPdf2, weight);
    if (std::abs(pdf->id1()) == 21)
      pdf_g->Fill(logPdf1, weight);
    if (std::abs(pdf->id2()) == 21)
      pdf_g->Fill(logPdf2, weight);
  }

  //Looping through the VERTICES in the event
  HepMC::GenEvent::vertex_const_iterator vrtxBegin = myGenEvent->vertices_begin();
  HepMC::GenEvent::vertex_const_iterator vrtxEnd = myGenEvent->vertices_end();
  unsigned int nvtx(0);
  for (HepMC::GenEvent::vertex_const_iterator vrtxIt = vrtxBegin; vrtxIt != vrtxEnd; ++vrtxIt) {
    HepMC::GenVertex const *vrtx = *vrtxIt;
    outVrtxPtclNumber->Fill(vrtx->particles_out_size(), weight);
    //std::cout << "all " << vrtx->particles_out_size() << '\n';

    if (nvtx == 0) {
      vrtxZ->Fill(vrtx->point3d().z(), weight);
      vrtxRadius->Fill(vrtx->point3d().perp(), weight);
    }
    HepMC::GenVertex::particles_out_const_iterator vrtxPtclBegin = vrtx->particles_out_const_begin();
    HepMC::GenVertex::particles_out_const_iterator vrtxPtclEnd = vrtx->particles_out_const_end();
    outVrtxStablePtclNum = 0;
    for (HepMC::GenVertex::particles_out_const_iterator vrtxPtclIt = vrtxPtclBegin; vrtxPtclIt != vrtxPtclEnd;
         ++vrtxPtclIt) {
      HepMC::GenParticle const *vrtxPtcl = *vrtxPtclIt;
      if (vrtxPtcl->status() == 1) {
        ++outVrtxStablePtclNum;
        //std::cout << "stable " << outVrtxStablePtclNum << '\n';
      }
    }
    outVrtxStablePtclNumber->Fill(outVrtxStablePtclNum, weight);
    nvtx++;
  }  //vertices

  HepMC::GenEvent::particle_const_iterator ptclBegin = myGenEvent->particles_begin();
  HepMC::GenEvent::particle_const_iterator ptclEnd = myGenEvent->particles_end();
  for (HepMC::GenEvent::particle_const_iterator ptclIt = ptclBegin; ptclIt != ptclEnd; ++ptclIt) {
    HepMC::GenParticle const *ptcl = *ptclIt;
    int Id = ptcl->pdg_id();  // std::cout << Id << '\n';
    float Log_p = log10(ptcl->momentum().rho());
    double charge = 999.;  // for the charge it's needed a HepPDT method
    int status = ptcl->status();
    const HepPDT::ParticleData *PData = fPDGTable->particle(HepPDT::ParticleID(Id));
    if (PData == nullptr) {
      //        std::cout << "Unknown id = " << Id << '\n';
      ++unknownPDTNum;
    } else
      charge = PData->charge();

    genPtclStatus->Fill((float)status, weight);

    if (ptcl->status() == 1) {
      ++stablePtclNum;
      stablePtclPhi->Fill(ptcl->momentum().phi() / CLHEP::degree,
                          weight);  //std::cout << ptcl->polarization().phi() << '\n';
      stablePtclEta->Fill(ptcl->momentum().pseudoRapidity(), weight);
      stablePtclCharge->Fill(charge, weight);  // std::cout << ptclData.charge() << '\n';
      stablePtclp->Fill(Log_p, weight);
      stablePtclpT->Fill(log10(ptcl->momentum().perp()), weight);
      if (charge != 0. && charge != 999.)
        ++stableChaNum;
      if (std::abs(Id) == 1)
        status1ShortLived->Fill(1, weight);
      if (std::abs(Id) == 2)
        status1ShortLived->Fill(2, weight);
      if (std::abs(Id) == 3)
        status1ShortLived->Fill(3, weight);
      if (std::abs(Id) == 4)
        status1ShortLived->Fill(4, weight);
      if (std::abs(Id) == 5)
        status1ShortLived->Fill(5, weight);
      if (std::abs(Id) == 6)
        status1ShortLived->Fill(6, weight);
      if (Id == 21)
        status1ShortLived->Fill(7, weight);
      if (std::abs(Id) == 15)
        status1ShortLived->Fill(8, weight);
      if (Id == 23)
        status1ShortLived->Fill(9, weight);
      if (std::abs(Id) == 24)
        status1ShortLived->Fill(10, weight);
      if (std::abs(Id) == 7 || std::abs(Id) == 8 || std::abs(Id) == 17 || (std::abs(Id) >= 25 && std::abs(Id) <= 99))
        status1ShortLived->Fill(11, weight);
      etotal += ptcl->momentum().e();
      pxtotal += ptcl->momentum().px();
      pytotal += ptcl->momentum().py();
      pztotal += ptcl->momentum().pz();
    }

    if (abs(Id) < 6 || abs(Id) == 22) {
      ++partonNum;
      partonpT->Fill(Log_p, weight);
    }

    bool indentified = false;
    for (unsigned int i = 0; i < particles.size(); i++) {
      if (particles.at(i).Fill(ptcl, weight)) {
        indentified = true;
        break;
      }
    }
    if (!indentified) {
      ++otherPtclNum;
      otherPtclMomentum->Fill(Log_p, weight);
    }
  }  //event particles

  // set a default sqrt(s) and then check in the event
  double ecms = 13000.;
  if (myGenEvent->valid_beam_particles()) {
    ecms = myGenEvent->beam_particles().first->momentum().e() + myGenEvent->beam_particles().second->momentum().e();
  }
  log10DeltaEcms->Fill(log10(fabs(etotal - ecms)), weight);
  DeltaEcms->Fill(etotal - ecms, weight);
  DeltaPx->Fill(pxtotal, weight);
  DeltaPy->Fill(pytotal, weight);
  DeltaPz->Fill(pztotal, weight);

  stablePtclNumber->Fill(log10(stablePtclNum + 0.1), weight);
  stableChaNumber->Fill(log10(stableChaNum + 0.1), weight);
  otherPtclNumber->Fill(log10(otherPtclNum + 0.1), weight);
  unknownPDTNumber->Fill(log10(unknownPDTNum + 0.1), weight);
  //
  partonNumber->Fill(partonNum, weight);
  for (unsigned int i = 0; i < particles.size(); i++) {
    particles.at(i).FillCount(weight);
  };

}  //analyze