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