g4simulation/g4jets/JetHepMCLoader.cc

0001 // $Id: $
0002
0003 /*!
0004  * \file JetHepMCLoader.C
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "JetHepMCLoader.h"
0012
0013 #include <jetbase/JetContainer.h>  // for JetContainer
0014 #include <jetbase/Jetv1.h>
0015
0016 #include <phhepmc/PHHepMCGenEvent.h>
0017 #include <phhepmc/PHHepMCGenEventMap.h>
0018
0019 #include <fun4all/Fun4AllBase.h>  // for Fun4AllBase::VERBOSITY_A_LOT
0020 #include <fun4all/Fun4AllHistoManager.h>
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/Fun4AllServer.h>
0023 #include <fun4all/SubsysReco.h>  // for SubsysReco
0024
0025 #include <phool/PHCompositeNode.h>
0026 #include <phool/PHIODataNode.h>
0027 #include <phool/PHNode.h>  // for PHNode
0028 #include <phool/PHNodeIterator.h>
0029 #include <phool/PHObject.h>  // for PHObject
0030 #include <phool/getClass.h>
0031 #include <phool/phool.h>  // for PHWHERE
0032
0033 #include <TAxis.h>  // for TAxis
0034 #include <TH1.h>
0035 #include <TH2.h>
0036 #include <TNamed.h>  // for TNamed
0037
0038 #include <HepMC/GenEvent.h>  // for GenEvent, GenEvent::particl...
0039 #include <HepMC/GenParticle.h>   // for GenParticle
0040 #include <HepMC/SimpleVector.h>  // for FourVector
0041 #include <HepMC/Units.h>         // for conversion_factor, GEV
0042
0043 #include <algorithm>  // for max
0044 #include <cassert>
0045 #include <iostream>
0046
0047 JetHepMCLoader::JetHepMCLoader(const std::string &jetInputCategory)
0048   : SubsysReco("JetHepMCLoader_" + jetInputCategory)
0049   , m_jetInputCategory(jetInputCategory)
0050
0051 {
0052   return;
0053 }
0054
0055 int JetHepMCLoader::InitRun(PHCompositeNode *topNode)
0056 {
0057   if (m_saveQAPlots)
0058   {
0059     Fun4AllHistoManager *hm = getHistoManager();
0060     assert(hm);
0061
0062     const int n_bins = 1 + m_jetSrc.size();
0063
0064     TH1D *h = new TH1D("hNormalization",  //
0065                        "Normalization;Items;Summed quantity", n_bins, .5, n_bins + .5);
0066     int i = 1;
0067     h->GetXaxis()->SetBinLabel(i++, "Event count");
0068     for (const hepmc_jet_src &src : m_jetSrc)
0069     {
0070       h->GetXaxis()->SetBinLabel(i++, (std::string("SubEvent ") + src.m_name).c_str());
0071     }
0072     h->GetXaxis()->LabelsOption("v");
0073     hm->registerHisto(h);
0074
0075     for (const hepmc_jet_src &src : m_jetSrc)
0076     {
0077       hm->registerHisto(
0078           new TH2F((std::string("hJetEtEta_") + src.m_name).c_str(),  //
0079                    ("Jet distribution from " + src.m_name + ";Jet #eta;Jet E_{T} [GeV]").c_str(),
0080                    40, -4., 4.,
0081                    100, 0, 100));
0082     }
0083   }  //   if (m_saveQAPlots)
0084
0085   return CreateNodes(topNode);
0086 }
0087
0088 int JetHepMCLoader::process_event(PHCompositeNode *topNode)
0089 {
0090   // For pile-up simulation: define GenEventMap
0091   PHHepMCGenEventMap *genevtmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0092
0093   if (!genevtmap)
0094   {
0095     static bool once = true;
0096
0097     if (once && Verbosity())
0098     {
0099       once = false;
0100
0101       std::cout << "HepMCNodeReader::process_event - No PHHepMCGenEventMap node. Do not perform HepMC->Geant4 input" << std::endl;
0102     }
0103
0104     return Fun4AllReturnCodes::DISCARDEVENT;
0105   }
0106
0107   if (m_saveQAPlots)
0108   {
0109     Fun4AllHistoManager *hm = getHistoManager();
0110     assert(hm);
0111     TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto("hNormalization"));
0112     assert(h_norm);
0113     h_norm->Fill("Event count", 1);
0114   }  //   if (m_saveQAPlots)
0115
0116   for (const hepmc_jet_src &src : m_jetSrc)
0117   {
0118     JetContainer *jets = findNode::getClass<JetContainer>(topNode, src.m_name);
0119     assert(jets);
0120
0121     jets->set_algo(src.m_algorithmID);
0122     jets->set_par(src.m_parameter);
0123     jets->insert_src(Jet::HEPMC_IMPORT);
0124
0125     PHHepMCGenEvent *genevt = genevtmap->get(src.m_embeddingID);
0126
0127     if (genevt == nullptr)
0128     {
0129       continue;
0130     }
0131
0132     HepMC::GenEvent *evt = genevt->getEvent();
0133     if (!evt)
0134     {
0135       std::cout << PHWHERE << " no evt pointer under HEPMC Node found:";
0136       genevt->identify();
0137       return Fun4AllReturnCodes::ABORTEVENT;
0138     }
0139
0140     assert(genevt->get_embedding_id() == src.m_embeddingID);
0141
0142     TH2F *hjet = nullptr;
0143
0144     if (m_saveQAPlots)
0145     {
0146       Fun4AllHistoManager *hm = getHistoManager();
0147       assert(hm);
0148       TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto("hNormalization"));
0149       assert(h_norm);
0150       h_norm->Fill((std::string("SubEvent ") + src.m_name).c_str(), 1);
0151
0152       hjet = dynamic_cast<TH2F *>(hm->getHisto(std::string("hJetEtEta_") + src.m_name));
0153       assert(hjet);
0154
0155     }  //   if (m_saveQAPlots)
0156
0157     const double mom_factor = HepMC::Units::conversion_factor(evt->momentum_unit(), HepMC::Units::GEV);
0158
0159     for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin();
0160          p != evt->particles_end(); ++p)
0161     {
0162       HepMC::GenParticle *part = (*p);
0163
0164       assert(part);
0165       if (Verbosity() >= VERBOSITY_A_LOT)
0166       {
0167         part->print();
0168       }
0169
0170       if (part->status() == src.m_tagStatus && part->pdg_id() == src.m_tagPID)
0171       {
0172         Jet *jet = jets->add_jet();  // returns a new Jetv2
0173
0174         jet->set_px(part->momentum().px() * mom_factor);
0175         jet->set_py(part->momentum().py() * mom_factor);
0176         jet->set_pz(part->momentum().pz() * mom_factor);
0177         jet->set_e(part->momentum().e() * mom_factor);
0178
0179         jet->insert_comp(Jet::HEPMC_IMPORT, part->barcode());
0180         if (hjet)
0181         {
0182           hjet->Fill(jet->get_eta(), jet->get_et());
0183         }
0184
0185       }  //       if (part->status() == src.m_tagStatus and part->pdg_id() == src.m_tagPID)
0186     }
0187   }  //  for (const hepmc_jet_src &src : m_jetSrc)
0188
0189   return Fun4AllReturnCodes::EVENT_OK;
0190 }
0191
0192 int JetHepMCLoader::End(PHCompositeNode * /*topNode*/)
0193 {
0194   if (m_saveQAPlots)
0195   {
0196     Fun4AllHistoManager *hm = getHistoManager();
0197     assert(hm);
0198
0199     std::cout << "JetHepMCLoader::End - saving QA histograms to " << Name() + ".root" << std::endl;
0200     hm->dumpHistos(Name() + ".root", "RECREATE");
0201   }
0202
0203   return Fun4AllReturnCodes::EVENT_OK;
0204 }
0205
0206 void JetHepMCLoader::addJet(
0207     const std::string &name,
0208     int embeddingID,
0209     Jet::ALGO algorithm,
0210     double parameter,
0211     int tagPID,
0212     int tagStatus)
0213 {
0214   std::string algorithmName = "Undefined_Jet_Algorithm";
0215
0216   switch (algorithm)
0217   {
0218   case Jet::ANTIKT:
0219     algorithmName = "ANTIKT";
0220     break;
0221
0222   case Jet::KT:
0223     algorithmName = "KT";
0224     break;
0225
0226   case Jet::CAMBRIDGE:
0227     algorithmName = "CAMBRIDGE";
0228     break;
0229
0230   default:
0231
0232     break;
0233   }
0234
0235   hepmc_jet_src src{name, embeddingID, algorithmName, algorithm, parameter, tagPID, tagStatus};
0236
0237   m_jetSrc.push_back(src);
0238 }
0239
0240 int JetHepMCLoader::CreateNodes(PHCompositeNode *topNode)
0241 {
0242   PHNodeIterator iter(topNode);
0243
0244   // Looking for the DST node
0245   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0246   if (!dstNode)
0247   {
0248     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0249     return Fun4AllReturnCodes::ABORTRUN;
0250   }
0251
0252   for (const hepmc_jet_src &src : m_jetSrc)
0253   {
0254     // Create the AntiKt node if required
0255     PHCompositeNode *AlgoNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", src.m_algorithmName));
0256     if (!AlgoNode)
0257     {
0258       AlgoNode = new PHCompositeNode(src.m_algorithmName);
0259       dstNode->addNode(AlgoNode);
0260     }
0261
0262     // Create the Input node if required
0263     PHCompositeNode *InputNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", m_jetInputCategory));
0264     if (!InputNode)
0265     {
0266       InputNode = new PHCompositeNode(m_jetInputCategory);
0267       AlgoNode->addNode(InputNode);
0268     }
0269
0270     JetContainer *jets = findNode::getClass<JetContainer>(topNode, src.m_name);
0271     if (!jets)
0272     {
0273       jets = new JetContainer();
0274       PHIODataNode<PHObject> *JetContainerNode = new PHIODataNode<PHObject>(jets, src.m_name, "PHObject");
0275       InputNode->addNode(JetContainerNode);
0276     }
0277   }
0278
0279   return Fun4AllReturnCodes::EVENT_OK;
0280 }
0281
0282 Fun4AllHistoManager *
0283 JetHepMCLoader::getHistoManager()
0284 {
0285   std::string histname(Name() + "_HISTOS");
0286
0287   Fun4AllServer *se = Fun4AllServer::instance();
0288   Fun4AllHistoManager *hm = se->getHistoManager(histname);
0289
0290   if (! hm)
0291   {
0292     std::cout
0293         << "TPCDataStreamEmulator::get_HistoManager - Making Fun4AllHistoManager " << histname
0294         << std::endl;
0295     hm = new Fun4AllHistoManager(histname);
0296     se->registerHistoManager(hm);
0297   }
0298
0299   assert(hm);
0300
0301   return hm;
0302 }