analysis/HERWIG_UE/HerwigSTARUECheck.cc

0001 #include "HerwigSTARUECheck.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <phool/PHCompositeNode.h>
0005 #include <phool/getClass.h>
0006 #include <phool/phool.h>
0007
0008 #include <phhepmc/PHHepMCGenEvent.h>
0009 #include <phhepmc/PHHepMCGenEventMap.h>
0010
0011 #include <g4main/PHG4Particle.h>
0012 #include <g4main/PHG4TruthInfoContainer.h>
0013
0014 #include <HepMC/GenEvent.h>
0015 #include <HepMC/GenParticle.h>
0016 #include <HepMC/SimpleVector.h>
0017 #include <HepMC/Units.h>
0018
0019 #include <fastjet/ClusterSequence.hh>
0020 #include <fastjet/JetDefinition.hh>
0021 #include <fastjet/PseudoJet.hh>
0022
0023 #include <TDatabasePDG.h>
0024 #include <TFile.h>
0025 #include <TH1D.h>
0026 #include <TH2D.h>
0027 #include <TMath.h>
0028 #include <TParticlePDG.h>
0029
0030 #include <cmath>
0031 #include <cstdlib>
0032 #include <iostream>
0033 #include <limits>
0034 #include <vector>
0035
0036 namespace
0037 {
0038   // STAR PRD 101, 052004 leading-jet pT binning (HEPData d01-x01-y03)
0039   const int NBINS = 8;
0040   const double BINS[NBINS + 1] = {5., 7., 9., 11., 15., 20., 25., 35., 45.};
0041
0042   const int NBINS_R04 = 9;
0043   const double BINS_R04[NBINS_R04 + 1] = {5., 7., 9., 11., 15., 21., 26., 32.5, 40.5, 63.5};
0044
0045   // STAR PLB 637 (2006) pion spectra binning
0046   const int NBINS_PI = 23;
0047   const double BINS_PI[NBINS_PI + 1] = {0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10.0};
0048 }  // namespace
0049
0050 // ====================================================================
0051 //  UEHistSet -- one full set of observables per particle source
0052 // ====================================================================
0053
0054 TH2 *HerwigSTARUECheck::UEHistSet::mk2(const std::string &name, const std::string &title,
0055                                        int nx, const double *xbins,
0056                                        int ny, double ylo, double yhi)
0057 {
0058   TH2D *h = new TH2D(name.c_str(), title.c_str(), nx, xbins, ny, ylo, yhi);
0059   h->Sumw2();
0060   all.push_back(h);
0061   return h;
0062 }
0063
0064 TH1 *HerwigSTARUECheck::UEHistSet::mk1(const std::string &name, const std::string &title,
0065                                        int nx, const double *xbins)
0066 {
0067   TH1D *h = new TH1D(name.c_str(), title.c_str(), nx, xbins);
0068   h->Sumw2();
0069   all.push_back(h);
0070   return h;
0071 }
0072
0073 TH1 *HerwigSTARUECheck::UEHistSet::mk1(const std::string &name, const std::string &title,
0074                                        int nx, double xlo, double xhi)
0075 {
0076   TH1D *h = new TH1D(name.c_str(), title.c_str(), nx, xlo, xhi);
0077   h->Sumw2();
0078   all.push_back(h);
0079   return h;
0080 }
0081
0082 void HerwigSTARUECheck::UEHistSet::book(const std::string &tag, const std::string &label)
0083 {
0084   const std::string t = "_" + tag;          // name suffix, e.g. "_hepmc"
0085   const std::string l = " (" + label + ")";  // title suffix, e.g. " (HepMC)"
0086
0087   dens_trans_02 = mk2("p_dens_trans_pt02" + t,
0088                       "Transverse #LTdN_{ch}/d#etad#phi#GT, p_{T}>0.2" + l +
0089                           ";leading jet p_{T} [GeV/c];#LTdN_{ch}/d#etad#phi#GT",
0090                       NBINS, BINS, 101, -0.5, 100.5);
0091   dens_toward_02 = mk2("p_dens_toward_pt02" + t,
0092                        "Toward #LTdN_{ch}/d#etad#phi#GT, p_{T}>0.2" + l +
0093                            ";leading jet p_{T} [GeV/c];#LTdN_{ch}/d#etad#phi#GT",
0094                        NBINS, BINS, 101, -0.5, 100.5);
0095   dens_away_02 = mk2("p_dens_away_pt02" + t,
0096                      "Away #LTdN_{ch}/d#etad#phi#GT, p_{T}>0.2" + l +
0097                          ";leading jet p_{T} [GeV/c];#LTdN_{ch}/d#etad#phi#GT",
0098                      NBINS, BINS, 101, -0.5, 100.5);
0099   dens_trans_05 = mk2("p_dens_trans_pt05" + t,
0100                       "Transverse #LTdN_{ch}/d#etad#phi#GT, p_{T}>0.5" + l +
0101                           ";leading jet p_{T} [GeV/c];#LTdN_{ch}/d#etad#phi#GT",
0102                       NBINS, BINS, 101, -0.5, 100.5);
0103
0104   avgpt_trans_02 = mk2("p_avgpt_trans_pt02" + t,
0105                        "Transverse #LTp_{T,ch}#GT, p_{T}>0.2" + l +
0106                            ";leading jet p_{T} [GeV/c];#LTp_{T}#GT [GeV/c]",
0107                        NBINS, BINS, 1000, 0.0, 70.0);
0108   avgpt_toward_02 = mk2("p_avgpt_toward_pt02" + t,
0109                         "Toward #LTp_{T,ch}#GT, p_{T}>0.2" + l +
0110                             ";leading jet p_{T} [GeV/c];#LTp_{T}#GT [GeV/c]",
0111                         NBINS, BINS, 1000, 0.0, 70.0);
0112   avgpt_away_02 = mk2("p_avgpt_away_pt02" + t,
0113                       "Away #LTp_{T,ch}#GT, p_{T}>0.2" + l +
0114                           ";leading jet p_{T} [GeV/c];#LTp_{T}#GT [GeV/c]",
0115                       NBINS, BINS, 1000, 0.0, 70.0);
0116   avgpt_trans_05 = mk2("p_avgpt_trans_pt05" + t,
0117                        "Transverse #LTp_{T,ch}#GT, p_{T}>0.5" + l +
0118                            ";leading jet p_{T} [GeV/c];#LTp_{T}#GT [GeV/c]",
0119                        NBINS, BINS, 1000, 0.0, 70.0);
0120
0121   sumET_trans_02 = mk2("p_sumET_trans_02" + t,
0122                        "Transverse #LT#Sigma E_{T}#GT, p_{T}>0.2" + l +
0123                            ";leading jet p_{T} [GeV/c];#LT#Sigma E_{T}/#delta#eta#delta#phi#GT [GeV]",
0124                        NBINS, BINS, 1000, 0.0, 70.0);
0125   sumET_toward_02 = mk2("p_sumET_toward_02" + t,
0126                         "Toward #LT#Sigma E_{T}#GT, p_{T}>0.2" + l +
0127                             ";leading jet p_{T} [GeV/c];#LT#Sigma E_{T}/#delta#eta#delta#phi#GT [GeV]",
0128                         NBINS, BINS, 1000, 0.0, 70.0);
0129   sumET_away_02 = mk2("p_sumET_away_02" + t,
0130                       "Away #LT#Sigma E_{T}#GT, p_{T}>0.2" + l +
0131                           ";leading jet p_{T} [GeV/c];#LT#Sigma E_{T}/#delta#eta#delta#phi#GT [GeV]",
0132                       NBINS, BINS, 1000, 0.0, 70.0);
0133   sumET_trans_05 = mk2("p_sumET_trans_05" + t,
0134                        "Transverse #LT#Sigma E_{T}#GT, p_{T}>0.5" + l +
0135                            ";leading jet p_{T} [GeV/c];#LT#Sigma E_{T}/#delta#eta#delta#phi#GT [GeV]",
0136                        NBINS, BINS, 1000, 0.0, 70.0);
0137
0138   sumET_trans_R04 = mk2("p_sumET_trans_R04" + t,
0139                         "Transverse #LT#Sigma E_{T}/#delta#eta#delta#phi#GT, p_{T,ch}>0.5, p_{T,neut}>0.2" + l +
0140                             ";leading R=0.4 jet p_{T} [GeV/c];#LT#Sigma E_{T}/#delta#eta#delta#phi#GT [GeV]",
0141                         NBINS_R04, BINS_R04, 1000, 0.0, 70.0);
0142
0143   pi_spec = mk1("pip_spec" + t,
0144                 "#pi^{+} spectrum" + l +
0145                     ";p_{T} [GeV/c];#frac{1}{2#pi p_{T}} d^{2}N/dp_{T}dy [GeV^{-2}]",
0146                 NBINS_PI, BINS_PI);
0147
0148   leadjet_pt = mk1("h_leadjet_pt" + t,
0149                    "accepted leading jet p_{T}" + l + ";p_{T} [GeV/c];events",
0150                    80, 0., 80.);
0151   leadjet_pt_R04 = mk1("h_leadjet_pt_R04" + t,
0152                        "accepted leading R=0.4 jet p_{T}" + l + ";p_{T} [GeV/c];events",
0153                        80, 0., 80.);
0154   leadjet_eta = mk1("h_leadjet_eta" + t,
0155                     "accepted leading jet #eta" + l + ";#eta;events",
0156                     40, -1., 1.);
0157   dphi = mk1("h_dphi" + t,
0158              "|#Delta#phi(part, leading jet)|, p_{T}>0.2" + l + ";|#Delta#phi|;particles",
0159              60, 0., M_PI);
0160
0161   cutflow = mk1("h_cutflow" + t, "cutflow" + l + ";;events", 3, 0.5, 3.5);
0162   cutflow->GetXaxis()->SetBinLabel(1, "particles found");
0163   cutflow->GetXaxis()->SetBinLabel(2, "leading R=0.6 jet accepted");
0164   cutflow->GetXaxis()->SetBinLabel(3, "leading R=0.4 jet accepted");
0165 }
0166
0167 void HerwigSTARUECheck::UEHistSet::write()
0168 {
0169   for (TH1 *h : all)
0170   {
0171     if (h) h->Write();
0172   }
0173 }
0174
0175 // ====================================================================
0176 //  module
0177 // ====================================================================
0178
0179 HerwigSTARUECheck::HerwigSTARUECheck(const std::string &name,
0180                                      const std::string &outfile)
0181   : SubsysReco(name)
0182   , m_outfileName(outfile)
0183 {
0184 }
0185
0186 int HerwigSTARUECheck::Init(PHCompositeNode * /*topNode*/)
0187 {
0188   if (!m_doHepMC && !m_doG4Truth)
0189   {
0190     std::cerr << Name() << ": both particle sources are disabled, nothing to do" << std::endl;
0191     return Fun4AllReturnCodes::ABORTRUN;
0192   }
0193
0194   m_outfile = new TFile(m_outfileName.c_str(), "RECREATE");
0195
0196   if (m_doHepMC) m_hHepMC.book("hepmc", "HepMC");
0197   if (m_doG4Truth) m_hG4.book("g4", "G4 primary");
0198
0199   m_h_nevents = new TH1D("h_nevents", "cutflow;;events", 3, 0.5, 3.5);
0200   m_h_nevents->GetXaxis()->SetBinLabel(1, "processed");
0201   m_h_nevents->GetXaxis()->SetBinLabel(2, "HepMC found");
0202   m_h_nevents->GetXaxis()->SetBinLabel(3, "G4TruthInfo found");
0203
0204   const bool isHerwig = (m_sampleName.find("Herwig") != std::string::npos);
0205   for (int s = 0; s < 8; s++)
0206   {
0207     if (!isHerwig && m_sampleName == sampleNames[s])
0208     {
0209       sampleNumber = s;
0210       break;
0211     }
0212     if (isHerwig && s < 7 && m_sampleName == HerwigsampleNames[s])
0213     {
0214       sampleNumber = s;
0215       break;
0216     }
0217   }
0218
0219   if (sampleNumber == -999)
0220   {
0221     std::cerr << "sample number is still default. Sample name " << m_sampleName
0222               << " is not valid. Exiting" << std::endl;
0223     return Fun4AllReturnCodes::ABORTRUN;
0224   }
0225
0226   // one weight per event, shared by both branches so the comparison is
0227   // apples-to-apples
0228   m_weight = m_useCSWeights ? (isHerwig ? HerwigCS[sampleNumber] : cs[sampleNumber]) : 1.0;
0229
0230   return Fun4AllReturnCodes::EVENT_OK;
0231 }
0232
0233 // --------------------------------------------------------------------
0234 //  particle collection
0235 // --------------------------------------------------------------------
0236
0237 bool HerwigSTARUECheck::collect_hepmc(PHCompositeNode *topNode,
0238                                       std::vector<fastjet::PseudoJet> &parts,
0239                                       std::vector<int> &pids)
0240 {
0241   PHHepMCGenEventMap *genmap =
0242       findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0243   if (!genmap)
0244   {
0245     static bool once = true;
0246     if (once)
0247     {
0248       std::cout << Name() << ": no PHHepMCGenEventMap node found -- "
0249                 << "is this a generator-level DST?" << std::endl;
0250       once = false;
0251     }
0252     return false;
0253   }
0254
0255   PHHepMCGenEvent *genevt = genmap->get(m_embeddingId);
0256   if (!genevt && !genmap->empty())
0257   {
0258     genevt = genmap->begin()->second;  // fall back to first embedded event
0259   }
0260   if (!genevt || !genevt->getEvent()) return false;
0261
0262   HepMC::GenEvent *evt = genevt->getEvent();
0263
0264   // sPHENIX convention is GeV, but guard against MeV productions
0265   const double mom2gev = (evt->momentum_unit() == HepMC::Units::MEV) ? 1.e-3 : 1.0;
0266
0267   for (HepMC::GenEvent::particle_const_iterator it = evt->particles_begin();
0268        it != evt->particles_end(); ++it)
0269   {
0270     const HepMC::GenParticle *p = *it;
0271     if (p->status() != 1) continue;
0272
0273     const int pid = p->pdg_id();
0274     const int apid = std::abs(pid);
0275     if (!m_includeNeutrinos && (apid == 12 || apid == 14 || apid == 16)) continue;
0276
0277     const HepMC::FourVector &mom = p->momentum();
0278     fastjet::PseudoJet pj(mom.px() * mom2gev, mom.py() * mom2gev,
0279                           mom.pz() * mom2gev, mom.e() * mom2gev);
0280
0281     if (pj.pt() <= m_constPtMin) continue;
0282     if (std::fabs(pj.pseudorapidity()) >= m_partEtaMax) continue;
0283
0284     parts.push_back(pj);
0285     pids.push_back(pid);
0286   }
0287
0288   return true;
0289 }
0290
0291 bool HerwigSTARUECheck::collect_g4truth(PHCompositeNode *topNode,
0292                                         std::vector<fastjet::PseudoJet> &parts,
0293                                         std::vector<int> &pids)
0294 {
0295   PHG4TruthInfoContainer *truthinfo =
0296       findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0297   if (!truthinfo)
0298   {
0299     static bool once = true;
0300     if (once)
0301     {
0302       std::cout << PHWHERE << "PHG4TruthInfoContainer node is missing, "
0303                 << "can't collect G4 truth particles" << std::endl;
0304       once = false;
0305     }
0306     return false;
0307   }
0308
0309   PHG4TruthInfoContainer::Range range = truthinfo->GetSPHENIXPrimaryParticleRange();
0310   for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0311   {
0312     const PHG4Particle *truth = iter->second;
0313
0314     const int pid = truth->get_pid();
0315     const int apid = std::abs(pid);
0316     if (!m_includeNeutrinos && (apid == 12 || apid == 14 || apid == 16)) continue;
0317
0318     fastjet::PseudoJet pj(truth->get_px(), truth->get_py(), truth->get_pz(), truth->get_e());
0319
0320     if (pj.pt() <= m_constPtMin) continue;
0321     if (std::fabs(pj.pseudorapidity()) >= m_partEtaMax) continue;
0322
0323     parts.push_back(pj);
0324     pids.push_back(pid);
0325   }
0326
0327   return true;
0328 }
0329
0330 // --------------------------------------------------------------------
0331 //  analysis, run once per particle source
0332 // --------------------------------------------------------------------
0333
0334 void HerwigSTARUECheck::analyze(UEHistSet &H,
0335                                 const std::vector<fastjet::PseudoJet> &parts,
0336                                 const std::vector<int> &pids,
0337                                 double w)
0338 {
0339   // pion spectrum: every event, no jet requirement (as before)
0340   for (int p = 0; p < (int) parts.size(); p++)
0341   {
0342     if (pids[p] == 211 && std::fabs(parts[p].rap()) < 0.5)
0343     {
0344       H.pi_spec->Fill(parts[p].pt(), w / (2.0 * TMath::Pi() * parts[p].pt()));
0345     }
0346   }
0347
0348   if (parts.empty()) return;
0349   H.cutflow->Fill(1);
0350
0351   // ------------------------------------------------------------------
0352   // jet finding: anti-kT R = m_jetR (0.6) and R = 0.4 on ALL selected
0353   // final-state particles
0354   // ------------------------------------------------------------------
0355   fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, m_jetR);
0356   fastjet::ClusterSequence cseq(parts, jetdef);
0357   std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cseq.inclusive_jets(m_jetPtMin));
0358
0359   fastjet::JetDefinition jetdefR04(fastjet::antikt_algorithm, 0.4);
0360   fastjet::ClusterSequence cseqR04(parts, jetdefR04);
0361   std::vector<fastjet::PseudoJet> jetsR04 = fastjet::sorted_by_pt(cseqR04.inclusive_jets(m_jetPtMin));
0362
0363   const fastjet::PseudoJet *lead = nullptr;
0364   for (const auto &j : jets)
0365   {
0366     if (std::fabs(j.pseudorapidity()) < m_jetEtaMax && j.pt() > m_jetPtMin && j.pt() < m_jetPtMax)
0367     {
0368       lead = &j;
0369       break;  // jets are pT-sorted, so this is the hardest jet in acceptance
0370     }
0371   }
0372
0373   // pThat-slice stitching window
0374   const double loR06 = truthJet_min_pT[sampleNumber];
0375   const double hiR06 = (sampleNumber + 1 < 8) ? truthJet_min_pT[sampleNumber + 1]
0376                                               : std::numeric_limits<double>::max();
0377   const double loR04 = truthJet_min_pT_R04[sampleNumber];
0378   const double hiR04 = (sampleNumber + 1 < 8) ? truthJet_min_pT_R04[sampleNumber + 1]
0379                                               : std::numeric_limits<double>::max();
0380
0381   bool goodR06 = (lead && lead->pt() >= loR06 && lead->pt() < hiR06);
0382
0383   const fastjet::PseudoJet *leadR04 = nullptr;
0384   for (const auto &j : jetsR04)
0385   {
0386     if (std::fabs(j.pseudorapidity()) < 0.7 && j.pt() > m_jetPtMin && j.pt() < m_jetPtMax)
0387     {
0388       leadR04 = &j;
0389       break;
0390     }
0391   }
0392
0393   bool goodR04 = (leadR04 && leadR04->pt() >= loR04 && leadR04->pt() < hiR04);
0394
0395   if (!goodR06 && !goodR04) return;
0396
0397   if (goodR06)
0398   {
0399     H.cutflow->Fill(2);
0400     H.leadjet_pt->Fill(lead->pt(), w);
0401     H.leadjet_eta->Fill(lead->pseudorapidity(), w);
0402   }
0403   if (goodR04)
0404   {
0405     H.cutflow->Fill(3);
0406     H.leadjet_pt_R04->Fill(leadR04->pt(), w);
0407   }
0408
0409   // ------------------------------------------------------------------
0410   // region assignment and densities
0411   // ------------------------------------------------------------------
0412   int nToward = 0;
0413   int nAway = 0;
0414   int nTrans = 0;
0415   int nTrans05 = 0;
0416
0417   double sumET_Toward = 0.0;
0418   double sumET_Away = 0.0;
0419   double sumET_Trans = 0.0;
0420   double sumET_Trans05 = 0.0;
0421   double sumET_TransR04 = 0.0;
0422
0423   if (goodR06)
0424   {
0425     for (int p = 0; p < (int) parts.size(); p++)
0426     {
0427       const fastjet::PseudoJet &ch = parts[p];
0428       const double adphi = std::fabs(ch.delta_phi_to(*lead));
0429       H.dphi->Fill(adphi, w);
0430
0431       const bool charged = (pdg_charge(pids[p]) != 0.);
0432
0433       if (adphi < M_PI / 3.)  // Toward
0434       {
0435         if (charged) ++nToward;
0436         sumET_Toward += ch.Et();
0437         H.avgpt_toward_02->Fill(lead->pt(), ch.pt(), w);
0438       }
0439       else if (adphi > 2. * M_PI / 3.)  // Away
0440       {
0441         if (charged) ++nAway;
0442         sumET_Away += ch.Et();
0443         H.avgpt_away_02->Fill(lead->pt(), ch.pt(), w);
0444       }
0445       else  // Transverse: pi/3 < |dphi| < 2pi/3
0446       {
0447         if (charged) ++nTrans;
0448         sumET_Trans += ch.Et();
0449         H.avgpt_trans_02->Fill(lead->pt(), ch.pt(), w);
0450         if (ch.pt() > 0.5)
0451         {
0452           if (charged) ++nTrans05;
0453           sumET_Trans05 += ch.Et();
0454           H.avgpt_trans_05->Fill(lead->pt(), ch.pt(), w);
0455         }
0456       }
0457     }
0458   }
0459
0460   if (goodR04)
0461   {
0462     for (int p = 0; p < (int) parts.size(); p++)
0463     {
0464       const fastjet::PseudoJet &part = parts[p];
0465       const double adphiR04 = std::fabs(part.delta_phi_to(*leadR04));
0466       if (adphiR04 >= M_PI / 3. && adphiR04 <= 2. * M_PI / 3.)
0467       {
0468         const bool em = (pids[p] == 22 || std::abs(pids[p]) == 11);
0469         if ((!em && part.pt() > 0.5) || (em && part.pt() > 0.2)) sumET_TransR04 += part.Et();
0470       }
0471     }
0472   }
0473
0474   // area normalization: deta * dphi = (2*etaMax) * (2pi/3), identical for
0475   // all three regions (each spans a total of 2pi/3 in azimuth)
0476   const double area = (2. * m_partEtaMax) * (2. * M_PI / 3.);
0477   const double areaR04 = (2. * 1.1) * (2. * M_PI / 3.);
0478
0479   if (goodR06)
0480   {
0481     H.dens_toward_02->Fill(lead->pt(), nToward, w);
0482     H.dens_away_02->Fill(lead->pt(), nAway, w);
0483     H.dens_trans_02->Fill(lead->pt(), nTrans, w);
0484     H.dens_trans_05->Fill(lead->pt(), nTrans05, w);
0485
0486     H.sumET_toward_02->Fill(lead->pt(), sumET_Toward / area, w);
0487     H.sumET_away_02->Fill(lead->pt(), sumET_Away / area, w);
0488     H.sumET_trans_02->Fill(lead->pt(), sumET_Trans / area, w);
0489     H.sumET_trans_05->Fill(lead->pt(), sumET_Trans05 / area, w);
0490   }
0491   if (goodR04) H.sumET_trans_R04->Fill(leadR04->pt(), sumET_TransR04 / areaR04, w);
0492 }
0493
0494 // --------------------------------------------------------------------
0495
0496 int HerwigSTARUECheck::process_event(PHCompositeNode *topNode)
0497 {
0498   m_h_nevents->Fill(1);
0499
0500   if (m_doHepMC)
0501   {
0502     std::vector<fastjet::PseudoJet> parts;
0503     std::vector<int> pids;
0504     if (collect_hepmc(topNode, parts, pids))
0505     {
0506       m_h_nevents->Fill(2);
0507       analyze(m_hHepMC, parts, pids, m_weight);
0508     }
0509     else if (m_requireBoth)
0510     {
0511       return Fun4AllReturnCodes::ABORTEVENT;
0512     }
0513   }
0514
0515   if (m_doG4Truth)
0516   {
0517     std::vector<fastjet::PseudoJet> parts;
0518     std::vector<int> pids;
0519     if (collect_g4truth(topNode, parts, pids))
0520     {
0521       m_h_nevents->Fill(3);
0522       analyze(m_hG4, parts, pids, m_weight);
0523     }
0524     else if (m_requireBoth)
0525     {
0526       return Fun4AllReturnCodes::ABORTEVENT;
0527     }
0528   }
0529
0530   return Fun4AllReturnCodes::EVENT_OK;
0531 }
0532
0533 int HerwigSTARUECheck::End(PHCompositeNode * /*topNode*/)
0534 {
0535   m_outfile->cd();
0536
0537   if (m_doHepMC) m_hHepMC.write();
0538   if (m_doG4Truth) m_hG4.write();
0539   m_h_nevents->Write();
0540
0541   m_outfile->Close();
0542
0543   std::cout << Name() << ": wrote " << m_outfileName << " ("
0544             << (long long) m_h_nevents->GetBinContent(1) << " events processed";
0545   if (m_doHepMC)
0546   {
0547     std::cout << "; HepMC: " << (long long) m_hHepMC.cutflow->GetBinContent(2)
0548               << " with an accepted leading R=0.6 jet";
0549   }
0550   if (m_doG4Truth)
0551   {
0552     std::cout << "; G4 primary: " << (long long) m_hG4.cutflow->GetBinContent(2)
0553               << " with an accepted leading R=0.6 jet";
0554   }
0555   std::cout << ")" << std::endl;
0556
0557   return Fun4AllReturnCodes::EVENT_OK;
0558 }
0559
0560 double HerwigSTARUECheck::pdg_charge(int pid)
0561 {
0562   auto it = m_chargeCache.find(pid);
0563   if (it != m_chargeCache.end()) return it->second;
0564
0565   double q = 0.;
0566   TParticlePDG *pdgp = TDatabasePDG::Instance()->GetParticle(pid);
0567   if (pdgp)
0568   {
0569     q = pdgp->Charge() / 3.;  // TParticlePDG stores charge in |e|/3
0570   }
0571   else if (Verbosity() > 0)
0572   {
0573     std::cout << Name() << ": PDG id " << pid
0574               << " unknown to TDatabasePDG, treating as neutral" << std::endl;
0575   }
0576   m_chargeCache[pid] = q;
0577   return q;
0578 }