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