Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:21:24 

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