g4simulation/g4eval/PHG4DSTReader.cc

0001 // $Id: PHG4DSTReader.cc,v 1.11 2015/01/06 02:52:07 jinhuang Exp $
0002
0003 /*!
0004  * \file PHG4DSTReader.cc
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision: 1.11 $
0008  * \date $Date: 2015/01/06 02:52:07 $
0009  */
0010
0011 #include "PHG4DSTReader.h"
0012
0013 #include <g4main/PHG4HitContainer.h>
0014 #include <g4main/PHG4TruthInfoContainer.h>
0015
0016 #include <calobase/RawTower.h>
0017 #include <calobase/RawTowerContainer.h>
0018 #include <calobase/RawTowerv1.h>
0019
0020 #include <g4main/PHG4Hit.h>
0021 #include <g4main/PHG4HitEval.h>
0022 #include <g4main/PHG4Particle.h>
0023 #include <g4main/PHG4Particlev2.h>
0024 #include <g4main/PHG4VtxPoint.h>
0025 #include <g4main/PHG4VtxPointv1.h>
0026
0027 #include <jetbase/Jet.h>
0028 #include <jetbase/JetContainer.h>
0029 #include <jetbase/Jetv2.h>
0030
0031 #include <fun4all/Fun4AllReturnCodes.h>
0032 #include <fun4all/PHTFileServer.h>
0033 #include <fun4all/SubsysReco.h>
0034
0035 #include <phool/getClass.h>
0036
0037 #include <TClass.h>
0038 #include <TClonesArray.h>
0039 #include <TObject.h>
0040 #include <TTree.h>
0041
0042 #include <cassert>
0043 #include <format>
0044 #include <iostream>
0045 #include <set>
0046 #include <sstream>
0047 #include <utility>
0048
0049 using part_type = PHG4Particlev2;
0050 using hit_type = PHG4HitEval;
0051 using vertex_type = PHG4VtxPointv1;
0052 using RawTower_type = RawTowerv1;
0053 using PHPyJet_type = Jetv2;
0054
0055 PHG4DSTReader::PHG4DSTReader(const std::string &filename)
0056   : SubsysReco("PHG4DSTReader")
0057   , _out_file_name(filename)
0058 {
0059   // TODO Auto-generated constructor stub
0060 }
0061
0062 PHG4DSTReader::~PHG4DSTReader()
0063 {
0064   if (Verbosity() > 1)
0065   {
0066     std::cout << "PHG4DSTReader::destructor - Clean ups" << std::endl;
0067   }
0068   if (_T)
0069   {
0070     _T->ResetBranchAddresses();
0071   }
0072
0073   _records.clear();
0074 }
0075
0076 int PHG4DSTReader::Init(PHCompositeNode * /*unused*/)
0077 {
0078   const int arr_size = 100;
0079
0080   //  BOOST_FOREACH(std::string nodenam, _node_postfix)
0081   for (const auto &nodenam : _node_postfix)
0082   {
0083     const char *class_name = hit_type::Class()->GetName();
0084
0085     std::string hname = std::format("G4HIT_{}", nodenam);
0086     //      _node_name.push_back(hname);
0087     if (Verbosity() > 0)
0088     {
0089       std::cout << "PHG4DSTReader::Init - saving hits from node: " << hname << " - "
0090                 << class_name << std::endl;
0091     }
0092
0093     record rec;
0094     rec._cnt = 0;
0095     rec._name = hname;
0096     rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
0097     rec._arr_ptr = rec._arr.get();
0098     rec._type = record::typ_hit;
0099
0100     _records.push_back(rec);
0101
0102     nblocks++;
0103   }
0104
0105   if (!_tower_postfix.empty() && Verbosity() > 0)
0106   {
0107     std::cout << "PHG4DSTReader::Init - zero suppression for calorimeter towers = "
0108               << _tower_zero_sup << " GeV" << std::endl;
0109   }
0110   for (const auto &nodenam : _tower_postfix)
0111   {
0112     const char *class_name = RawTower_type::Class()->GetName();
0113
0114     std::string hname = std::format("TOWER_{}", nodenam);
0115     //      _node_name.push_back(hname);
0116     if (Verbosity() > 0)
0117     {
0118       std::cout << "PHG4DSTReader::Init - saving raw tower info from node: " << hname
0119                 << " - " << class_name << std::endl;
0120     }
0121     record rec;
0122     rec._cnt = 0;
0123     rec._name = hname;
0124     rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
0125     rec._arr_ptr = rec._arr.get();
0126     rec._type = record::typ_tower;
0127
0128     _records.push_back(rec);
0129
0130     nblocks++;
0131   }
0132
0133   for (const auto &hname : _jet_postfix)
0134   {
0135     const char *class_name = PHPyJet_type::Class()->GetName();
0136
0137     if (Verbosity() > 0)
0138     {
0139       std::cout << "PHG4DSTReader::Init - saving jets from node: " << hname << " - "
0140                 << class_name << std::endl;
0141     }
0142     record rec;
0143     rec._cnt = 0;
0144     rec._name = hname;
0145     rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
0146     rec._arr_ptr = rec._arr.get();
0147     rec._type = record::typ_jets;
0148
0149     _records.push_back(rec);
0150
0151     nblocks++;
0152   }
0153
0154   if (_save_particle)
0155   {
0156     // save particles
0157
0158     const char *class_name = part_type::Class()->GetName();
0159
0160     if (Verbosity() > 0)
0161     {
0162       std::cout << "PHG4DSTReader::Init - saving Particles node:" << class_name
0163                 << std::endl;
0164     }
0165     record rec;
0166     rec._cnt = 0;
0167     rec._name = "PHG4Particle";
0168     rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
0169     rec._arr_ptr = rec._arr.get();
0170     rec._type = record::typ_part;
0171
0172     _records.push_back(rec);
0173
0174     nblocks++;
0175   }
0176
0177   if (_save_vertex)
0178   {
0179     // save particles
0180
0181     const char *class_name = vertex_type::Class()->GetName();
0182
0183     if (Verbosity() > 0)
0184     {
0185       std::cout << "PHG4DSTReader::Init - saving vertex for particles under node:"
0186                 << class_name << std::endl;
0187     }
0188     record rec;
0189     rec._cnt = 0;
0190     rec._name = "PHG4VtxPoint";
0191     rec._arr = boost::make_shared<TClonesArray>(class_name, arr_size);
0192     rec._arr_ptr = rec._arr.get();
0193     rec._type = record::typ_vertex;
0194
0195     _records.push_back(rec);
0196
0197     nblocks++;
0198   }
0199   if (Verbosity() > 0)
0200   {
0201     std::cout << "PHG4DSTReader::Init - requested " << nblocks << " nodes" << std::endl;
0202   }
0203   build_tree();
0204
0205   return 0;
0206 }
0207
0208 void PHG4DSTReader::build_tree()
0209 {
0210   if (Verbosity() > 0)
0211   {
0212     std::cout << "PHG4DSTReader::build_tree - output to " << _out_file_name << std::endl;
0213   }
0214   static const int BUFFER_SIZE = 32000;
0215
0216   // open TFile
0217   PHTFileServer::open(_out_file_name, "RECREATE");
0218
0219   _T = new TTree("T", "PHG4DSTReader");
0220
0221   nblocks = 0;
0222   for (auto &rec : _records)
0223   {
0224     if (Verbosity() > 0)
0225     {
0226       std::cout << "PHG4DSTReader::build_tree - Add " << rec._name << std::endl;
0227     }
0228     const std::string name_cnt = "n_" + rec._name;
0229     const std::string name_cnt_desc = name_cnt + "/I";
0230     _T->Branch(name_cnt.c_str(), &(rec._cnt), name_cnt_desc.c_str(),
0231                BUFFER_SIZE);
0232     _T->Branch(rec._name.c_str(), &(rec._arr_ptr), BUFFER_SIZE, 99);
0233
0234     nblocks++;
0235   }
0236
0237   if (Verbosity() > 0)
0238   {
0239     std::cout << "PHG4DSTReader::build_tree - added " << nblocks << " nodes" << std::endl;
0240   }
0241   _T->SetAutoSave(16000);
0242 }
0243
0244 int PHG4DSTReader::process_event(PHCompositeNode *topNode)
0245 {
0246   _event++;
0247
0248   // clean ups
0249   _particle_set.clear();
0250   _vertex_set.clear();
0251
0252   PHG4TruthInfoContainer *truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0253   if (!truthInfoList)
0254   {
0255     if (_event < 2)
0256     {
0257       std::cout
0258           << "PHG4DSTReader::process_event - Error - can not find node G4TruthInfo. Quit processing!"
0259           << std::endl;
0260     }
0261     return Fun4AllReturnCodes::ABORTEVENT;
0262   }
0263
0264   for (auto &rec : _records)
0265   {
0266     rec._cnt = 0;
0267     assert(rec._arr.get() == rec._arr_ptr);
0268     assert(rec._arr.get());
0269     rec._arr->Clear();
0270
0271     if (rec._type == record::typ_hit)
0272     {
0273       if (Verbosity() >= 2)
0274       {
0275         std::cout << "PHG4DSTReader::process_event - processing " << rec._name
0276                   << std::endl;
0277       }
0278
0279       PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode,
0280                                                                     rec._name);
0281       if (!hits)
0282       {
0283         if (_event < 2)
0284         {
0285           std::cout
0286               << "PHG4DSTReader::process_event - Error - can not find node "
0287               << rec._name << std::endl;
0288         }
0289       }
0290       else
0291       {
0292         PHG4HitContainer::ConstRange hit_range = hits->getHits();
0293
0294         if (Verbosity() >= 2)
0295         {
0296           std::cout << "PHG4DSTReader::process_event - processing "
0297                     << rec._name << " and received " << hits->size() << " hits"
0298                     << std::endl;
0299         }
0300
0301         for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first;
0302              hit_iter != hit_range.second; hit_iter++)
0303         {
0304           PHG4Hit *hit = hit_iter->second;
0305
0306           //                  hit_type * hit = dynamic_cast<hit_type *>(hit_raw);
0307           //
0308           assert(hit);
0309
0310           new ((*(rec._arr.get()))[rec._cnt]) hit_type(hit);
0311
0312           hit_type *new_hit =
0313               dynamic_cast<hit_type *>(rec._arr.get()->At(rec._cnt));
0314           assert(new_hit);
0315
0316           //                  for (int i = 0; i < 2; i++)
0317           //                    {
0318           //                      new_hit->set_x(i, hit->get_x(i));
0319           //                      new_hit->set_y(i, hit->get_y(i));
0320           //                      new_hit->set_z(i, hit->get_z(i));
0321           //                      new_hit->set_t(i, hit->get_t(i));
0322           //                    }
0323           //                  new_hit->set_edep(hit->get_edep());
0324           //                  new_hit->set_layer(hit->get_layer());
0325           //                  new_hit->set_hit_id(hit->get_hit_id());
0326           //                  new_hit->set_trkid(hit->get_trkid());
0327
0328           //                  *new_hit = (*hit);
0329
0330           if (_load_active_particle)
0331           {
0332             _particle_set.insert(hit->get_trkid());
0333           }
0334
0335           if (Verbosity() >= 2)
0336           {
0337             std::cout << "PHG4DSTReader::process_event - processing "
0338                       << rec._name << " and hit " << hit->get_hit_id()
0339                       << " with track id " << hit->get_trkid() << std::endl;
0340           }
0341
0342           rec._cnt++;
0343         }
0344       }  // if (!hits)
0345     }    //      if (rec._type == record::typ_hit)
0346     else if (rec._type == record::typ_tower)
0347     {
0348       if (Verbosity() >= 2)
0349       {
0350         std::cout << "PHG4DSTReader::process_event - processing tower "
0351                   << rec._name << std::endl;
0352       }
0353
0354       RawTowerContainer *hits = findNode::getClass<RawTowerContainer>(
0355           topNode, rec._name);
0356       if (!hits)
0357       {
0358         if (_event < 2)
0359         {
0360           std::cout
0361               << "PHG4DSTReader::process_event - Error - can not find node "
0362               << rec._name << std::endl;
0363         }
0364       }
0365       else
0366       {
0367         RawTowerContainer::ConstRange hit_range = hits->getTowers();
0368
0369         if (Verbosity() >= 2)
0370         {
0371           std::cout << "PHG4DSTReader::process_event - processing "
0372                     << rec._name << " and received " << hits->size()
0373                     << " tower hits" << std::endl;
0374         }
0375
0376         for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
0377              hit_iter != hit_range.second; hit_iter++)
0378         {
0379           RawTower *hit_raw = hit_iter->second;
0380
0381           RawTower_type *hit = dynamic_cast<RawTower_type *>(hit_raw);
0382           //                  RawTower * hit = hit_iter->second;
0383
0384           assert(hit);
0385
0386           if (hit->get_energy() < _tower_zero_sup)
0387           {
0388             if (Verbosity() >= 2)
0389             {
0390               std::cout
0391                   << "PHG4DSTReader::process_event - suppress low energy tower hit "
0392                   << rec._name << " @ ("
0393                   //                            << hit->get_thetaMin()
0394                   //                            << ", " << hit->get_phiMin()
0395                   << "), Energy = " << hit->get_energy() << std::endl;
0396             }
0397
0398             continue;
0399           }
0400
0401           new ((*(rec._arr.get()))[rec._cnt]) RawTower_type();
0402
0403           if (Verbosity() >= 2)
0404           {
0405             std::cout
0406                 << "PHG4DSTReader::process_event - processing Tower hit "
0407                 << rec._name << " @ ("
0408                 //                        << hit->get_thetaMin() << ", "
0409                 //                        << hit->get_phiMin()
0410                 << "), Energy = " << hit->get_energy() << " - "
0411                 << rec._arr.get()->At(rec._cnt)->ClassName() << std::endl;
0412           }
0413
0414           //                  rec._arr->Print();
0415
0416           RawTower_type *new_hit =
0417               dynamic_cast<RawTower_type *>(rec._arr.get()->At(rec._cnt));
0418           assert(new_hit);
0419
0420           //                  //v1 copy
0421           //                  new_hit->get_bineta(hit->get_bineta());
0422           //                  new_hit->get_binphi(hit->get_binphi());
0423           //
0424           //                  //v2 copy
0425           //                  new_hit->set_thetaMin(hit->get_thetaMin());
0426           //                  new_hit->set_thetaSize(hit->get_thetaSize());
0427           //                  new_hit->set_phiMin(hit->get_phiMin());
0428           //                  new_hit->set_phiSize(hit->get_phiSize());
0429           //                  new_hit->set_zMin(hit->get_zMin());
0430           //                  new_hit->set_zSize(hit->get_zSize());
0431
0432           *new_hit = (*hit);
0433
0434           rec._cnt++;
0435         }
0436       }  // if (!hits)
0437     }    //      if (rec._type == record::typ_hit)
0438     else if (rec._type == record::typ_jets)
0439     {
0440       //          std::std::cout
0441       //              << "PHG4DSTReader::AddJet - Error - temp. disabled until jet added back to sPHENIX software"
0442       //              << std::std::endl;
0443       //
0444       if (Verbosity() >= 2)
0445       {
0446         std::cout << "PHG4DSTReader::process_event - processing jets "
0447                   << rec._name << std::endl;
0448       }
0449
0450       JetContainer *hits = findNode::getClass<JetContainer>(topNode, rec._name);
0451       if (!hits)
0452       {
0453         if (_event < 2)
0454         {
0455           std::cout
0456               << "PHG4DSTReader::process_event - Error - can not find node "
0457               << rec._name << std::endl;
0458         }
0459       }
0460       else
0461       {
0462         if (Verbosity() >= 2)
0463         {
0464           std::cout << "PHG4DSTReader::process_event - processing "
0465                     << rec._name << " and received " << hits->size() << " jets"
0466                     << std::endl;
0467         }
0468
0469         // for every recojet
0470         for (auto *hit_raw : *hits)
0471         {
0472           /* Jet *hit_raw = iter.second; */
0473
0474           if (Verbosity() >= 2)
0475           {
0476             std::cout << "PHG4DSTReader::process_event - processing jet "
0477                       << rec._name << " @ (" << hit_raw->get_eta() << ", "
0478                       << hit_raw->get_phi() << "), pT = " << hit_raw->get_pt()
0479                       << " - with raw type " << hit_raw->ClassName() << std::endl;
0480           }
0481
0482           PHPyJet_type *hit = dynamic_cast<PHPyJet_type *>(hit_raw);
0483
0484           assert(hit);
0485
0486           new ((*(rec._arr.get()))[rec._cnt]) PHPyJet_type();
0487
0488           PHPyJet_type *new_hit =
0489               dynamic_cast<PHPyJet_type *>(rec._arr.get()->At(rec._cnt));
0490           assert(new_hit);
0491
0492           *new_hit = (Jetv2) (*hit);
0493
0494           rec._cnt++;
0495         }
0496       }  // if (!hits)
0497     }    //      if (rec._type == record::typ_hit)
0498     else if (rec._type == record::typ_part)
0499     {
0500       if (_load_all_particle)
0501       {
0502         static bool once = true;
0503         if (once)
0504         {
0505           std::cout
0506               << "PHG4DSTReader::process_event - will load all particle from G4TruthInfo"
0507               << std::endl;
0508
0509           once = false;
0510         }
0511
0512         for (auto particle_iter : truthInfoList->GetMap())
0513         {
0514           _particle_set.insert(particle_iter.first);
0515         }
0516
0517       }  //          if (_load_all_particle)
0518       else
0519       {
0520         static bool once = true;
0521         if (once)
0522         {
0523           if (Verbosity() > 0)
0524           {
0525             std::cout << "PHG4DSTReader::process_event - will load primary particle from G4TruthInfo" << std::endl;
0526           }
0527           once = false;
0528         }
0529
0530         PHG4TruthInfoContainer::ConstRange primary_range =
0531             truthInfoList->GetPrimaryParticleRange();
0532
0533         for (PHG4TruthInfoContainer::ConstIterator particle_iter =
0534                  primary_range.first;
0535              particle_iter != primary_range.second;
0536              ++particle_iter)
0537         {
0538           _particle_set.insert(particle_iter->first);
0539
0540         }  //          if (_load_all_particle) else
0541       }
0542       for (int i : _particle_set)
0543       {
0544         auto *particle = truthInfoList->GetParticle(i);
0545         if (!particle)
0546         {
0547           std::cout
0548               << "PHG4DSTReader::process_event - ERROR - can not find particle/track ID "
0549               << i << " in G4TruthInfo" << std::endl;
0550
0551           continue;
0552         }
0553
0554         add_particle(rec, particle);
0555       }  // for(PartSet_t::const_iterator i = _particle_set.begin();i!=_particle_set.end();i++)
0556
0557     }  //      if (rec._type == record::typ_part)
0558     else if (rec._type == record::typ_vertex)
0559     {
0560       static bool once = true;
0561       if (once)
0562       {
0563         if (Verbosity() > 0)
0564         {
0565           std::cout << "PHG4DSTReader::process_event - will load vertex from G4TruthInfo" << std::endl;
0566         }
0567         once = false;
0568       }
0569
0570       for (int i : _vertex_set)
0571       {
0572         PHG4VtxPoint *v = truthInfoList->GetVtx(i);
0573         if (!v)
0574         {
0575           std::cout
0576               << "PHG4DSTReader::process_event - ERROR - can not find vertex ID "
0577               << i << " in G4TruthInfo" << std::endl;
0578
0579           continue;
0580         }
0581
0582         new ((*(rec._arr.get()))[rec._cnt]) vertex_type();
0583
0584         if (Verbosity() >= 2)
0585         {
0586           std::cout << "PHG4DSTReader::process_event - saving vertex id "
0587                     << i << std::endl;
0588         }
0589
0590         vertex_type *new_v =
0591             static_cast<vertex_type *>(rec._arr.get()->At(rec._cnt));
0592         assert(new_v);
0593
0594         new_v->set_x(v->get_x());
0595         new_v->set_y(v->get_y());
0596         new_v->set_z(v->get_z());
0597         new_v->set_t(v->get_t());
0598         new_v->set_id(v->get_id());
0599
0600         rec._cnt++;
0601       }  // else if (rec._type == record::typ_vertex)
0602
0603     }  //          if (_load_all_particle)
0604
0605   }  //  for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
0606
0607   if (_T)
0608   {
0609     _T->Fill();
0610   }
0611
0612   return 0;
0613 }  //  for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
0614
0615 void PHG4DSTReader::add_particle(PHG4DSTReader::record &rec, PHG4Particle *part)
0616 {
0617   assert(part);
0618
0619   //              if (Verbosity() >= 2)
0620   //                std::cout << "PHG4DSTReader::process_event - Particle type is "
0621   //                    << p_raw->GetName() << " - " << p_raw->ClassName()
0622   //                    << " get_track_id = " << p_raw->get_track_id() << std::endl;
0623
0624   //              part_type * part = dynamic_cast<part_type *>(p_raw);
0625
0626   //              assert(part);
0627
0628   new ((*(rec._arr.get()))[rec._cnt]) part_type();
0629
0630   part_type *new_part = static_cast<part_type *>(rec._arr.get()->At(rec._cnt));
0631   assert(new_part);
0632
0633   new_part->set_track_id(part->get_track_id());
0634   new_part->set_vtx_id(part->get_vtx_id());
0635   new_part->set_parent_id(part->get_parent_id());
0636   new_part->set_primary_id(part->get_primary_id());
0637   new_part->set_name(part->get_name());
0638   new_part->set_pid(part->get_pid());
0639   new_part->set_px(part->get_px());
0640   new_part->set_py(part->get_py());
0641   new_part->set_pz(part->get_pz());
0642   new_part->set_e(part->get_e());
0643
0644   _vertex_set.insert(part->get_vtx_id());
0645
0646   rec._cnt++;
0647 }
0648
0649 int PHG4DSTReader::End(PHCompositeNode * /*topNode*/)
0650 {
0651   if (Verbosity() > 1)
0652   {
0653     std::cout << "PHG4DSTReader::End - Clean ups" << std::endl;
0654   }
0655
0656   if (_T)
0657   {
0658     PHTFileServer::cd(_out_file_name);
0659     _T->Write();
0660     _T->ResetBranchAddresses();
0661   }
0662
0663   _records.clear();
0664
0665   return 0;
0666 }