Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:18:48 

0001 #include "QAG4SimulationTruthDecay.h"
0002 
0003 #include <phool/PHNode.h>
0004 #include <qautils/QAHistManagerDef.h>
0005 
0006 #include <fun4all/Fun4AllHistoManager.h>
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <fun4all/SubsysReco.h>
0009 
0010 #include <decayfinder/DecayFinder.h>
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/getClass.h>
0013 
0014 #include <CLHEP/Vector/LorentzVector.h>
0015 
0016 #include <TH1.h>
0017 #include <TString.h>
0018 
0019 #include <algorithm>
0020 #include <cassert>
0021 #include <cmath>
0022 #include <cstdlib>
0023 #include <iostream>
0024 #include <iterator>
0025 #include <ostream>
0026 #include <string>
0027 #include <vector>
0028 
0029 #include <cmath>
0030 
0031 static int candidateCounter = 0;
0032 
0033 /*
0034  *  QA module to check that decayers obey laws of physics
0035  *  Can also measure geometrical acceptance in eta and pT
0036  *  (Useful for cross-section measurements and similar)
0037  *  Authors: Cameron Dean and Thomas Marshall
0038  *  Date: July 2022
0039  */
0040 
0041 //____________________________________________________________________________..
0042 QAG4SimulationTruthDecay::QAG4SimulationTruthDecay(const std::string &name)
0043   : SubsysReco(name)
0044   , m_nTracks(0)
0045   , m_pt_min(0.2)
0046   , m_eta_min(-1.1)
0047   , m_eta_max(1.1)
0048   , m_decay_pdg_id(0)
0049   , m_truth_info(nullptr)
0050   , m_g4particle(nullptr)
0051   , m_df_module_name("myFinder")
0052   , m_outfile_name("outputData.root")
0053   , m_outfile(nullptr)
0054   , m_tree(nullptr)
0055   , m_write_nTuple(true)
0056   , m_write_QAHists(true)
0057 {
0058 }
0059 
0060 //____________________________________________________________________________..
0061 int QAG4SimulationTruthDecay::Init(PHCompositeNode *topNode)
0062 {
0063   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0064   assert(hm);
0065 
0066   TH1 *h(nullptr);
0067 
0068   int const m_mother_PDG_ID_nBins = (int) (m_mother_PDG_ID_max - m_mother_PDG_ID_min) * 1.2;
0069   int const m_daughter_PDG_ID_nBins = (int) (m_daughter_PDG_ID_max - m_daughter_PDG_ID_min) * 1.2;
0070 
0071   h = new TH1I(TString(get_histo_prefix()) + "mother_PDG_ID",  //
0072                ";Mother PDG ID;Entries", m_mother_PDG_ID_nBins, m_mother_PDG_ID_min, m_mother_PDG_ID_max);
0073 
0074   hm->registerHisto(h);
0075   h = new TH1F(TString(get_histo_prefix()) + "mother_mass",  //
0076                ";Mother Mass [GeV/c^{2}];Entries", m_mass_nBins, m_mass_min, m_mass_max);
0077   hm->registerHisto(h);
0078   h = new TH1F(TString(get_histo_prefix()) + "daughter_sum_mass",  //
0079                ";Daughter Sum Mass [GeV/c^{2}];Entries", m_mass_nBins, m_mass_min, m_mass_max);
0080   hm->registerHisto(h);
0081   h = new TH1F(TString(get_histo_prefix()) + "mother_decayLength",  //
0082                ";Mother Decay Length [cm];Entries", m_decayLength_nBins, m_decayLength_min, m_decayLength_max);
0083   hm->registerHisto(h);
0084   h = new TH1F(TString(get_histo_prefix()) + "mother_decayTime",  //
0085                ";Mother Decay Time [s];Entries", m_decayTime_nBins, m_decayTime_min, m_decayTime_max);
0086   hm->registerHisto(h);
0087   h = new TH1F(TString(get_histo_prefix()) + "mother_px",  //
0088                ";Mother p_{x} [GeV/c];Entries", 100, 0, 10);
0089   hm->registerHisto(h);
0090   h = new TH1F(TString(get_histo_prefix()) + "mother_py",  //
0091                ";Mother p_{y} [GeV/c];Entries", 100, 0, 10);
0092   hm->registerHisto(h);
0093   h = new TH1F(TString(get_histo_prefix()) + "mother_pz",  //
0094                ";Mother p_{z} [GeV/c];Entries", 100, 0, 10);
0095   hm->registerHisto(h);
0096   h = new TH1F(TString(get_histo_prefix()) + "mother_pE",  //
0097                ";Mother p_{E} [GeV];Entries", 100, 0, 10);
0098   hm->registerHisto(h);
0099   h = new TH1F(TString(get_histo_prefix()) + "mother_pT",  //
0100                ";Mother p_{T} [GeV/c];Entries", 100, 0, 10);
0101   hm->registerHisto(h);
0102   h = new TH1F(TString(get_histo_prefix()) + "mother_eta",  //
0103                ";Mother #eta;Entries", 1000, -10, 10);
0104   hm->registerHisto(h);
0105 
0106   for (unsigned int i = 0; i < 4; ++i)
0107   {
0108     std::string const track_number = "track_" + std::to_string(i + 1);
0109 
0110     h = new TH1I(TString(get_histo_prefix()) + TString(track_number) + "_PDG_ID",  //
0111                  ";Track PDG ID;Entries", m_daughter_PDG_ID_nBins, m_daughter_PDG_ID_min, m_daughter_PDG_ID_max);
0112     hm->registerHisto(h);
0113 
0114     /*
0115 
0116     h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_px",  //
0117                  ";Track p_{x} [GeV/c];Entries", 50, 0, 5);
0118     hm->registerHisto(h);
0119     h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_py",  //
0120                  ";Track p_{y} [GeV/c];Entries", 50, 0, 5);
0121     hm->registerHisto(h);
0122     h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_pz",  //
0123                  ";Track p_{z} [GeV/c];Entries", 50, 0, 5);
0124     hm->registerHisto(h);
0125     h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_pE",  //
0126                  ";Track p_{E} [GeV];Entries", 50, 0, 5);
0127     hm->registerHisto(h);
0128 
0129     */
0130 
0131     h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_pT",  //
0132                  ";Track pT [GeV/c];Entries", 50, 0, 5);
0133     hm->registerHisto(h);
0134     h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_eta",  //
0135                  ";Track Eta;Entries", 1000, -10, 10);
0136     hm->registerHisto(h);
0137     // h = new TH1F(TString(get_histo_prefix()) + TString(track_number) + "_mass",  //
0138     //              ";Track Mass [GeV/c^{2}];Entries", 100, 0, 3);
0139     // hm->registerHisto(h);
0140   }
0141 
0142   h = new TH1F(TString(get_histo_prefix()) + "delta_px",  //
0143                ";#delta p_{x} [GeV/c];Entries", 100, -3, 3);
0144   hm->registerHisto(h);
0145   h = new TH1F(TString(get_histo_prefix()) + "delta_py",  //
0146                ";#delta p_{y} [GeV/c];Entries", 100, -3, 3);
0147   hm->registerHisto(h);
0148   h = new TH1F(TString(get_histo_prefix()) + "delta_pz",  //
0149                ";#delta p_{z} [GeV/c];Entries", 100, -3, 3);
0150   hm->registerHisto(h);
0151   h = new TH1F(TString(get_histo_prefix()) + "delta_pE",  //
0152                ";#delta p_{E} [GeV];Entries", 100, -3, 3);
0153   hm->registerHisto(h);
0154 
0155   h = new TH1I(TString(get_histo_prefix()) + "accept_px_1percent",  //
0156                ";Accept p_{x} 1pcnt;Entries", 2, 0, 1);
0157   hm->registerHisto(h);
0158   h = new TH1I(TString(get_histo_prefix()) + "accept_py_1percent",  //
0159                ";Accept p_{y} 1pcnt;Entries", 2, 0, 1);
0160   hm->registerHisto(h);
0161   h = new TH1I(TString(get_histo_prefix()) + "accept_pz_1percent",  //
0162                ";Accept p_{z} 1pcnt;Entries", 2, 0, 1);
0163   hm->registerHisto(h);
0164   h = new TH1I(TString(get_histo_prefix()) + "accept_pE_1percent",  //
0165                ";Accept p_{E} 1pcnt;Entries", 2, 0, 1);
0166   hm->registerHisto(h);
0167 
0168   h = new TH1I(TString(get_histo_prefix()) + "accept_px_5percent",  //
0169                ";Accept p_{x} 5pcnt;Entries", 2, 0, 1);
0170   hm->registerHisto(h);
0171   h = new TH1I(TString(get_histo_prefix()) + "accept_py_5percent",  //
0172                ";Accept p_{y} 5pcnt;Entries", 2, 0, 1);
0173   hm->registerHisto(h);
0174   h = new TH1I(TString(get_histo_prefix()) + "accept_pz_5percent",  //
0175                ";Accept p_{z} 5pcnt;Entries", 2, 0, 1);
0176   hm->registerHisto(h);
0177   h = new TH1I(TString(get_histo_prefix()) + "accept_pE_5percent",  //
0178                ";Accept p_{E} 5pcnt;Entries", 2, 0, 1);
0179   hm->registerHisto(h);
0180 
0181   h = new TH1I(TString(get_histo_prefix()) + "accept_px_15percent",  //
0182                ";Accept p_{x} 15pcnt;Entries", 2, 0, 1);
0183   hm->registerHisto(h);
0184   h = new TH1I(TString(get_histo_prefix()) + "accept_py_15percent",  //
0185                ";Accept p_{y} 15pcnt;Entries", 2, 0, 1);
0186   hm->registerHisto(h);
0187   h = new TH1I(TString(get_histo_prefix()) + "accept_pz_15percent",  //
0188                ";Accept p_{z} 15pcnt;Entries", 2, 0, 1);
0189   hm->registerHisto(h);
0190   h = new TH1I(TString(get_histo_prefix()) + "accept_pE_15percent",  //
0191                ";Accept p_{E} 15pcnt;Entries", 2, 0, 1);
0192   hm->registerHisto(h);
0193 
0194   h = new TH1I(TString(get_histo_prefix()) + "accept_pT",  //
0195                ";Accept p_{T};Entries", 2, 0, 1);
0196   hm->registerHisto(h);
0197   h = new TH1I(TString(get_histo_prefix()) + "accept_eta",  //
0198                ";Accept #eta;Entries", 2, 0, 1);
0199   hm->registerHisto(h);
0200 
0201   assert(topNode);
0202   return Fun4AllReturnCodes::EVENT_OK;
0203 }
0204 
0205 //____________________________________________________________________________..
0206 int QAG4SimulationTruthDecay::process_event(PHCompositeNode *topNode)
0207 {
0208   resetValues();
0209   CLHEP::HepLorentzVector daughterSumLV;
0210 
0211   Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
0212   assert(hm);
0213 
0214   TH1I *h_mother_PDG_ID = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "mother_PDG_ID"));
0215   assert(h_mother_PDG_ID);
0216   TH1F *h_mother_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_mass"));
0217   assert(h_mother_mass);
0218   TH1F *h_daughter_sum_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "daughter_sum_mass"));
0219   assert(h_daughter_sum_mass);
0220   TH1F *h_mother_decayLength = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_decayLength"));
0221   assert(h_mother_decayLength);
0222   TH1F *h_mother_decayTime = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_decayTime"));
0223   assert(h_mother_decayTime);
0224   TH1F *h_mother_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_px"));
0225   assert(h_mother_px);
0226   TH1F *h_mother_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_py"));
0227   assert(h_mother_py);
0228   TH1F *h_mother_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_pz"));
0229   assert(h_mother_pz);
0230   TH1F *h_mother_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_pE"));
0231   assert(h_mother_pE);
0232   TH1F *h_mother_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_pT"));
0233   assert(h_mother_pT);
0234   TH1F *h_mother_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "mother_eta"));
0235   assert(h_mother_eta);
0236 
0237   TH1F *h_delta_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "delta_px"));
0238   assert(h_delta_px);
0239   TH1F *h_delta_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "delta_py"));
0240   assert(h_delta_py);
0241   TH1F *h_delta_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "delta_pz"));
0242   assert(h_delta_pz);
0243   TH1F *h_delta_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "delta_pE"));
0244   assert(h_delta_pE);
0245   TH1I *h_accept_px_1percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_px_1percent"));
0246   assert(h_accept_px_1percent);
0247   TH1I *h_accept_py_1percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_py_1percent"));
0248   assert(h_accept_py_1percent);
0249   TH1I *h_accept_pz_1percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_pz_1percent"));
0250   assert(h_accept_pz_1percent);
0251   TH1I *h_accept_pE_1percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_pE_1percent"));
0252   assert(h_accept_pE_1percent);
0253   TH1I *h_accept_px_5percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_px_5percent"));
0254   assert(h_accept_px_5percent);
0255   TH1I *h_accept_py_5percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_py_5percent"));
0256   assert(h_accept_py_5percent);
0257   TH1I *h_accept_pz_5percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_pz_5percent"));
0258   assert(h_accept_pz_5percent);
0259   TH1I *h_accept_pE_5percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_pE_5percent"));
0260   assert(h_accept_pE_5percent);
0261   TH1I *h_accept_px_15percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_px_15percent"));
0262   assert(h_accept_px_15percent);
0263   TH1I *h_accept_py_15percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_py_15percent"));
0264   assert(h_accept_py_15percent);
0265   TH1I *h_accept_pz_15percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_pz_15percent"));
0266   assert(h_accept_pz_15percent);
0267   TH1I *h_accept_pE_15percent = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_pE_15percent"));
0268   assert(h_accept_pE_15percent);
0269   TH1I *h_accept_pT = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_pT"));
0270   assert(h_accept_pT);
0271   TH1I *h_accept_eta = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "accept_eta"));
0272   assert(h_accept_eta);
0273 
0274   TH1I *h_track_1_PDG_ID = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "track_1_PDG_ID"));
0275   assert(h_track_1_PDG_ID);
0276   /*
0277   TH1F *h_track_1_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_px"));
0278   assert(h_track_1_px);
0279   TH1F *h_track_1_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_py"));
0280   assert(h_track_1_py);
0281   TH1F *h_track_1_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_pz"));
0282   assert(h_track_1_pz);
0283   TH1F *h_track_1_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_pE"));
0284   assert(h_track_1_pE);
0285   */
0286   TH1F *h_track_1_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_pT"));
0287   assert(h_track_1_pT);
0288   TH1F *h_track_1_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_eta"));
0289   assert(h_track_1_eta);
0290   // TH1F *h_track_1_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_1_mass"));
0291   // assert(h_track_1_mass);
0292 
0293   TH1I *h_track_2_PDG_ID = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "track_2_PDG_ID"));
0294   assert(h_track_2_PDG_ID);
0295   /*
0296   TH1F *h_track_2_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_px"));
0297   assert(h_track_2_px);
0298   TH1F *h_track_2_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_py"));
0299   assert(h_track_2_py);
0300   TH1F *h_track_2_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_pz"));
0301   assert(h_track_2_pz);
0302   TH1F *h_track_2_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_pE"));
0303   assert(h_track_2_pE);
0304   */
0305   TH1F *h_track_2_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_pT"));
0306   assert(h_track_2_pT);
0307   TH1F *h_track_2_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_eta"));
0308   assert(h_track_2_eta);
0309   // TH1F *h_track_2_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_2_mass"));
0310   // assert(h_track_2_mass);
0311 
0312   TH1I *h_track_3_PDG_ID = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "track_3_PDG_ID"));
0313   assert(h_track_3_PDG_ID);
0314   /*
0315   TH1F *h_track_3_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_px"));
0316   assert(h_track_3_px);
0317   TH1F *h_track_3_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_py"));
0318   assert(h_track_3_py);
0319   TH1F *h_track_3_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_pz"));
0320   assert(h_track_3_pz);
0321   TH1F *h_track_3_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_pE"));
0322   assert(h_track_3_pE);
0323   */
0324   TH1F *h_track_3_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_pT"));
0325   assert(h_track_3_pT);
0326   TH1F *h_track_3_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_eta"));
0327   assert(h_track_3_eta);
0328   // TH1F *h_track_3_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_3_mass"));
0329   // assert(h_track_3_mass);
0330 
0331   TH1I *h_track_4_PDG_ID = dynamic_cast<TH1I *>(hm->getHisto(get_histo_prefix() + "track_4_PDG_ID"));
0332   assert(h_track_4_PDG_ID);
0333   /*
0334   TH1F *h_track_4_px = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_px"));
0335   assert(h_track_4_px);
0336   TH1F *h_track_4_py = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_py"));
0337   assert(h_track_4_py);
0338   TH1F *h_track_4_pz = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_pz"));
0339   assert(h_track_4_pz);
0340   TH1F *h_track_4_pE = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_pE"));
0341   assert(h_track_4_pE);
0342   */
0343   TH1F *h_track_4_pT = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_pT"));
0344   assert(h_track_4_pT);
0345   TH1F *h_track_4_eta = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_eta"));
0346   assert(h_track_4_eta);
0347   // TH1F *h_track_4_mass = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "track_4_mass"));
0348   // assert(h_track_4_mass);
0349 
0350   ++m_event_number;
0351   if (m_decay_pdg_id == 0)
0352   {
0353     getMotherPDG(topNode);
0354   }
0355   std::vector<int> motherBarcodes = getDecayFinderMothers(topNode);
0356 
0357   if (motherBarcodes.size() == 1)
0358   {
0359     m_truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0360     if (!m_truth_info)
0361     {
0362       std::cout << "QAG4SimulationTruthDecay: Missing node G4TruthInfo" << std::endl;
0363       return Fun4AllReturnCodes::ABORTEVENT;
0364     }
0365 
0366     unsigned int trackCounter = 0;
0367     float mother_x = 0;
0368     float mother_y = 0;
0369     float mother_z = 0;
0370     float daughter_x = 0;
0371     float daughter_y = 0;
0372     float daughter_z = 0;
0373 
0374     PHG4TruthInfoContainer::ConstRange const range = m_truth_info->GetParticleRange();
0375     for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0376          iter != range.second; ++iter)
0377     {
0378       m_g4particle = iter->second;
0379       if (std::find(std::begin(motherBarcodes), std::end(motherBarcodes),
0380                     m_g4particle->get_barcode()) != std::end(motherBarcodes) &&
0381           abs(m_g4particle->get_pid()) == abs(m_decay_pdg_id))
0382       {
0383         m_mother_pdg_id = m_g4particle->get_pid();
0384         m_mother_barcode = m_g4particle->get_barcode();
0385         m_mother_px = m_g4particle->get_px();
0386         m_mother_py = m_g4particle->get_py();
0387         m_mother_pz = m_g4particle->get_pz();
0388         m_mother_pE = m_g4particle->get_e();
0389         CLHEP::HepLorentzVector const motherLV(m_mother_px, m_mother_py, m_mother_pz, m_mother_pE);
0390         m_mother_pT = motherLV.perp();
0391         m_mother_eta = motherLV.pseudoRapidity();
0392         m_mother_mass = motherLV.m();
0393 
0394         m_delta_px += m_mother_px;
0395         m_delta_py += m_mother_py;
0396         m_delta_pz += m_mother_pz;
0397         m_delta_pE += m_mother_pE;
0398 
0399         PHG4VtxPoint *mother_vtx = m_truth_info->GetVtx(m_g4particle->get_vtx_id());
0400         mother_x = mother_vtx->get_x();
0401         mother_y = mother_vtx->get_y();
0402         mother_z = mother_vtx->get_z();
0403 
0404         ++candidateCounter;
0405       }
0406 
0407       if (m_g4particle->get_parent_id() != 0)
0408       {
0409         PHG4Particle *mother = m_truth_info->GetParticle(m_g4particle->get_parent_id());
0410         if (std::find(std::begin(motherBarcodes), std::end(motherBarcodes),
0411                       mother->get_barcode()) != std::end(motherBarcodes) &&
0412             abs(mother->get_pid()) == abs(m_decay_pdg_id))
0413         {
0414           m_track_pdg_id[trackCounter] = m_g4particle->get_pid();
0415           m_track_mother_barcode[trackCounter] = mother->get_barcode();
0416           m_track_px[trackCounter] = m_g4particle->get_px();
0417           m_track_py[trackCounter] = m_g4particle->get_py();
0418           m_track_pz[trackCounter] = m_g4particle->get_pz();
0419           m_track_pE[trackCounter] = m_g4particle->get_e();
0420           CLHEP::HepLorentzVector const daughterLV(m_track_px[trackCounter], m_track_py[trackCounter], m_track_pz[trackCounter], m_track_pE[trackCounter]);
0421           daughterSumLV += daughterLV;
0422           m_track_pT[trackCounter] = daughterLV.perp();
0423           m_track_eta[trackCounter] = daughterLV.pseudoRapidity();
0424           m_track_mass[trackCounter] = daughterLV.m();
0425 
0426           m_delta_px -= m_track_px[trackCounter];
0427           m_delta_py -= m_track_py[trackCounter];
0428           m_delta_pz -= m_track_pz[trackCounter];
0429           m_delta_pE -= m_track_pE[trackCounter];
0430 
0431           PHG4VtxPoint *daughter_vtx = m_truth_info->GetVtx(m_g4particle->get_vtx_id());
0432           daughter_x = daughter_vtx->get_x();
0433           daughter_y = daughter_vtx->get_y();
0434           daughter_z = daughter_vtx->get_z();
0435 
0436           if (m_track_pT[trackCounter] < m_pt_min)
0437           {
0438             m_accept_pT = false;
0439           }
0440           bool const in_eta_range = isInRange(m_track_eta[trackCounter], m_eta_min, m_eta_max);
0441           if (!in_eta_range)
0442           {
0443             m_accept_eta = false;
0444           }
0445 
0446           ++trackCounter;
0447         }
0448       }
0449     }
0450 
0451     m_daughter_sum_mass = daughterSumLV.m();
0452 
0453     float const diff_percent_px = std::fabs(m_delta_px / m_mother_px) * 100.;
0454     float const diff_percent_py = std::fabs(m_delta_py / m_mother_py) * 100.;
0455     float const diff_percent_pz = std::fabs(m_delta_pz / m_mother_pz) * 100.;
0456     float const diff_percent_pE = std::fabs(m_delta_pE / m_mother_pE) * 100.;
0457 
0458     m_accept_px_1percent = diff_percent_px <= 1. ? true : false;
0459     m_accept_py_1percent = diff_percent_py <= 1. ? true : false;
0460     m_accept_pz_1percent = diff_percent_pz <= 1. ? true : false;
0461     m_accept_pE_1percent = diff_percent_pE <= 1. ? true : false;
0462 
0463     m_accept_px_5percent = diff_percent_px <= 5. ? true : false;
0464     m_accept_py_5percent = diff_percent_py <= 5. ? true : false;
0465     m_accept_pz_5percent = diff_percent_pz <= 5. ? true : false;
0466     m_accept_pE_5percent = diff_percent_pE <= 5. ? true : false;
0467 
0468     m_accept_px_15percent = diff_percent_px <= 15. ? true : false;
0469     m_accept_py_15percent = diff_percent_py <= 15. ? true : false;
0470     m_accept_pz_15percent = diff_percent_pz <= 15. ? true : false;
0471     m_accept_pE_15percent = diff_percent_pE <= 15. ? true : false;
0472 
0473     m_mother_decayLength = sqrt(pow(daughter_x - mother_x, 2) + pow(daughter_y - mother_y, 2) + pow(daughter_z - mother_z, 2));
0474     float const mother_p = sqrt(pow(m_mother_px, 2) + pow(m_mother_py, 2) + pow(m_mother_pz, 2));
0475     m_mother_decayTime = m_mother_mass * m_mother_decayLength / mother_p;
0476 
0477     if (m_write_nTuple)
0478     {
0479       if (candidateCounter == 1)
0480       {
0481         initializeBranches();
0482       }
0483       m_tree->Fill();
0484     }
0485     if (m_write_QAHists)
0486     {
0487       h_mother_PDG_ID->Fill(m_mother_pdg_id);
0488       h_mother_mass->Fill(m_mother_mass);
0489       h_daughter_sum_mass->Fill(m_daughter_sum_mass);
0490       h_mother_decayLength->Fill(m_mother_decayLength);
0491       h_mother_decayTime->Fill(m_mother_decayTime);
0492       h_mother_px->Fill(m_mother_px);
0493       h_mother_py->Fill(m_mother_py);
0494       h_mother_pz->Fill(m_mother_pz);
0495       h_mother_pE->Fill(m_mother_pE);
0496       h_mother_pT->Fill(m_mother_pT);
0497       h_mother_eta->Fill(m_mother_eta);
0498       h_delta_px->Fill(m_delta_px);
0499       h_delta_py->Fill(m_delta_py);
0500       h_delta_pz->Fill(m_delta_pz);
0501       h_delta_pE->Fill(m_delta_pE);
0502       h_accept_px_1percent->Fill(m_accept_px_1percent);
0503       h_accept_py_1percent->Fill(m_accept_py_1percent);
0504       h_accept_pz_1percent->Fill(m_accept_pz_1percent);
0505       h_accept_pE_1percent->Fill(m_accept_pE_1percent);
0506       h_accept_px_5percent->Fill(m_accept_px_5percent);
0507       h_accept_py_5percent->Fill(m_accept_py_5percent);
0508       h_accept_pz_5percent->Fill(m_accept_pz_5percent);
0509       h_accept_pE_5percent->Fill(m_accept_pE_5percent);
0510       h_accept_px_15percent->Fill(m_accept_px_15percent);
0511       h_accept_py_15percent->Fill(m_accept_py_15percent);
0512       h_accept_pz_15percent->Fill(m_accept_pz_15percent);
0513       h_accept_pE_15percent->Fill(m_accept_pE_15percent);
0514       h_accept_pT->Fill(m_accept_pT);
0515       h_accept_eta->Fill(m_accept_eta);
0516       if (m_nTracks >= 2)
0517       {
0518         h_track_1_PDG_ID->Fill(m_track_pdg_id[0]);
0519         /*
0520         h_track_1_px->Fill(m_track_px[0]);
0521         h_track_1_py->Fill(m_track_py[0]);
0522         h_track_1_pz->Fill(m_track_pz[0]);
0523         h_track_1_pE->Fill(m_track_pE[0]);
0524         */
0525         h_track_1_pT->Fill(m_track_pT[0]);
0526         h_track_1_eta->Fill(m_track_eta[0]);
0527         // h_track_1_mass->Fill(m_track_mass[0]);
0528         h_track_2_PDG_ID->Fill(m_track_pdg_id[1]);
0529         /*
0530         h_track_2_px->Fill(m_track_px[1]);
0531         h_track_2_py->Fill(m_track_py[1]);
0532         h_track_2_pz->Fill(m_track_pz[1]);
0533         h_track_2_pE->Fill(m_track_pE[1]);
0534         */
0535         h_track_2_pT->Fill(m_track_pT[1]);
0536         h_track_2_eta->Fill(m_track_eta[1]);
0537         // h_track_2_mass->Fill(m_track_mass[1]);
0538       }
0539       if (m_nTracks >= 3)
0540       {
0541         h_track_3_PDG_ID->Fill(m_track_pdg_id[2]);
0542         /*
0543         h_track_3_px->Fill(m_track_px[2]);
0544         h_track_3_py->Fill(m_track_py[2]);
0545         h_track_3_pz->Fill(m_track_pz[2]);
0546         h_track_3_pE->Fill(m_track_pE[2]);
0547         */
0548         h_track_3_pT->Fill(m_track_pT[2]);
0549         h_track_3_eta->Fill(m_track_eta[2]);
0550         // h_track_3_mass->Fill(m_track_mass[2]);
0551       }
0552       if (m_nTracks == 4)
0553       {
0554         h_track_4_PDG_ID->Fill(m_track_pdg_id[3]);
0555         /*
0556         h_track_4_px->Fill(m_track_px[3]);
0557         h_track_4_py->Fill(m_track_py[3]);
0558         h_track_4_pz->Fill(m_track_pz[3]);
0559         h_track_4_pE->Fill(m_track_pE[3]);
0560         */
0561         h_track_4_pT->Fill(m_track_pT[3]);
0562         h_track_4_eta->Fill(m_track_eta[3]);
0563         // h_track_4_mass->Fill(m_track_mass[3]);
0564       }
0565     }
0566   }
0567   else
0568   {
0569     std::cout << "You have more than one good decay in this event, this processing is not yet supported" << std::endl;
0570     return Fun4AllReturnCodes::ABORTEVENT;
0571   }
0572 
0573   return Fun4AllReturnCodes::EVENT_OK;
0574 }
0575 
0576 //____________________________________________________________________________..
0577 int QAG4SimulationTruthDecay::End(PHCompositeNode *topNode)
0578 {
0579   if (m_write_nTuple)
0580   {
0581     m_outfile->Write();
0582     m_outfile->Close();
0583     delete m_outfile;
0584   }
0585 
0586   assert(topNode);
0587 
0588   return Fun4AllReturnCodes::EVENT_OK;
0589 }
0590 
0591 void QAG4SimulationTruthDecay::initializeBranches()
0592 {
0593   m_outfile = new TFile(m_outfile_name.c_str(), "RECREATE");
0594   delete m_tree;
0595   m_tree = new TTree("QAG4SimulationTruthDecay", "QAG4SimulationTruthDecay");
0596   m_tree->OptimizeBaskets();
0597   m_tree->SetAutoSave(-5e6);  // Save the output file every 5MB
0598 
0599   m_tree->Branch("EventNumber", &m_event_number, "EventNumber/i");
0600   m_tree->Branch("mother_PDG_ID", &m_mother_pdg_id, "mother_PDG_ID/I");
0601   m_tree->Branch("mother_mass", &m_mother_mass, "mother_mass/F");
0602   m_tree->Branch("daughter_sum_mass", &m_daughter_sum_mass, "daughter_sum_mass/F");
0603   m_tree->Branch("mother_decayLength", &m_mother_decayLength, "mother_decayLength/F");
0604   m_tree->Branch("mother_decayTime", &m_mother_decayTime, "mother_decayTime/F");
0605   m_tree->Branch("mother_px", &m_mother_px, "mother_px/F");
0606   m_tree->Branch("mother_py", &m_mother_py, "mother_py/F");
0607   m_tree->Branch("mother_pz", &m_mother_pz, "mother_pz/F");
0608   m_tree->Branch("mother_pE", &m_mother_pE, "mother_pE/F");
0609   m_tree->Branch("mother_pT", &m_mother_pT, "mother_pT/F");
0610   m_tree->Branch("mother_eta", &m_mother_eta, "mother_eta/F");
0611   m_tree->Branch("mother_barcode", &m_mother_barcode, "mother_barcode/I");
0612 
0613   for (unsigned int i = 0; i < m_nTracks; ++i)
0614   {
0615     std::string const track_number = "track_" + std::to_string(i + 1);
0616 
0617     m_tree->Branch(TString(track_number) + "_PDG_ID", &m_track_pdg_id[i], TString(track_number) + "_PDG_ID/I");
0618     m_tree->Branch(TString(track_number) + "_px", &m_track_px[i], TString(track_number) + "_px/F");
0619     m_tree->Branch(TString(track_number) + "_py", &m_track_py[i], TString(track_number) + "_py/F");
0620     m_tree->Branch(TString(track_number) + "_pz", &m_track_pz[i], TString(track_number) + "_pz/F");
0621     m_tree->Branch(TString(track_number) + "_pE", &m_track_pE[i], TString(track_number) + "_pE/F");
0622     m_tree->Branch(TString(track_number) + "_pT", &m_track_pT[i], TString(track_number) + "_pT/F");
0623     m_tree->Branch(TString(track_number) + "_eta", &m_track_eta[i], TString(track_number) + "_eta/F");
0624     m_tree->Branch(TString(track_number) + "_mass", &m_track_mass[i], TString(track_number) + "_mass/F");
0625     m_tree->Branch(TString(track_number) + "_mother_barcode", &m_track_mother_barcode[i], TString(track_number) + "_mother_barcode/I");
0626   }
0627 
0628   m_tree->Branch("delta_px", &m_delta_px, "delta_px/F");
0629   m_tree->Branch("delta_py", &m_delta_py, "delta_py/F");
0630   m_tree->Branch("delta_pz", &m_delta_pz, "delta_pz/F");
0631   m_tree->Branch("delta_pE", &m_delta_pE, "delta_pE/F");
0632 
0633   m_tree->Branch("accept_px_1percent", &m_accept_px_1percent, "accept_px_1percent/O");
0634   m_tree->Branch("accept_py_1percent", &m_accept_py_1percent, "accept_py_1percent/O");
0635   m_tree->Branch("accept_pz_1percent", &m_accept_pz_1percent, "accept_pz_1percent/O");
0636   m_tree->Branch("accept_pE_1percent", &m_accept_pE_1percent, "accept_pE_1percent/O");
0637 
0638   m_tree->Branch("accept_px_5percent", &m_accept_px_5percent, "accept_px_5percent/O");
0639   m_tree->Branch("accept_py_5percent", &m_accept_py_5percent, "accept_py_5percent/O");
0640   m_tree->Branch("accept_pz_5percent", &m_accept_pz_5percent, "accept_pz_5percent/O");
0641   m_tree->Branch("accept_pE_5percent", &m_accept_pE_5percent, "accept_pE_5percent/O");
0642 
0643   m_tree->Branch("accept_px_15percent", &m_accept_px_15percent, "accept_px_15percent/O");
0644   m_tree->Branch("accept_py_15percent", &m_accept_py_15percent, "accept_py_15percent/O");
0645   m_tree->Branch("accept_pz_15percent", &m_accept_pz_15percent, "accept_pz_15percent/O");
0646   m_tree->Branch("accept_pE_15percent", &m_accept_pE_15percent, "accept_pE_15percent/O");
0647 
0648   m_tree->Branch("accept_pT", &m_accept_pT, "accept_pT/O");
0649   m_tree->Branch("accept_eta", &m_accept_eta, "accept_eta/O");
0650 }
0651 
0652 void QAG4SimulationTruthDecay::getMotherPDG(PHCompositeNode *topNode)
0653 {
0654   PHNodeIterator nodeIter(topNode);
0655 
0656   std::string const node_name = m_df_module_name + "_DecayMap";
0657 
0658   PHNode *findNode = dynamic_cast<PHNode *>(nodeIter.findFirst(node_name.c_str()));
0659   if (findNode)
0660   {
0661     m_decayMap = findNode::getClass<DecayFinderContainer_v1>(topNode, node_name.c_str());
0662   }
0663   else
0664   {
0665   }
0666 
0667   Decay decay = m_decayMap->begin()->second;
0668   m_decay_pdg_id = decay[0].second;
0669 }
0670 
0671 std::vector<int> QAG4SimulationTruthDecay::getDecayFinderMothers(PHCompositeNode *topNode)
0672 {
0673   std::vector<int> m_motherBarcodes;
0674 
0675   PHNodeIterator const nodeIter(topNode);
0676 
0677   std::string const node_name = m_df_module_name + "_DecayMap";
0678 
0679   m_decayMap = findNode::getClass<DecayFinderContainer_v1>(topNode, node_name.c_str());
0680 
0681   for (auto &iter : *m_decayMap)
0682   {
0683     Decay decay = iter.second;
0684     m_nTracks = decay.size() - 1;
0685     m_motherBarcodes.push_back(decay[0].first.second);
0686   }
0687 
0688   return m_motherBarcodes;
0689 }
0690 
0691 bool QAG4SimulationTruthDecay::isInRange(float min, float value, float max)
0692 {
0693   return min <= value && value <= max;
0694 }
0695 
0696 void QAG4SimulationTruthDecay::resetValues()
0697 {
0698   m_delta_px = 0;
0699   m_delta_py = 0;
0700   m_delta_pz = 0;
0701   m_delta_pE = 0;
0702 
0703   m_accept_eta = true;
0704   m_accept_pT = true;
0705 }
0706 
0707 std::string QAG4SimulationTruthDecay::get_histo_prefix()
0708 {
0709   return std::string("h_") + Name() + std::string("_") + std::string("_");
0710 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!