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 <TString.h>
0041 #include <TTree.h>
0042
0043 #include <cassert>
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 = Form("G4HIT_%s", nodenam.c_str());
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.size() && 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 = Form("TOWER_%s", nodenam.c_str());
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::get().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   //  const double significand = _event / TMath::Power(10, (int) (log10(_event)));
0247   //
0248   //  if (fmod(significand, 1.0) == 0 && significand <= 10)
0249   //    std::cout << "PHG4DSTReader::process_event - " << _event << std::endl;
0250   _event++;
0251
0252   // clean ups
0253   _particle_set.clear();
0254   _vertex_set.clear();
0255
0256   PHG4TruthInfoContainer *truthInfoList = findNode::getClass<
0257       PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0258   if (!truthInfoList)
0259   {
0260     if (_event < 2)
0261     {
0262       std::cout
0263           << "PHG4DSTReader::process_event - Error - can not find node G4TruthInfo. Quit processing!"
0264           << std::endl;
0265     }
0266     return Fun4AllReturnCodes::ABORTEVENT;
0267   }
0268
0269   for (auto &rec : _records)
0270   {
0271     rec._cnt = 0;
0272     assert(rec._arr.get() == rec._arr_ptr);
0273     assert(rec._arr.get());
0274     rec._arr->Clear();
0275
0276     if (rec._type == record::typ_hit)
0277     {
0278       if (Verbosity() >= 2)
0279       {
0280         std::cout << "PHG4DSTReader::process_event - processing " << rec._name
0281                   << std::endl;
0282       }
0283
0284       PHG4HitContainer *hits = findNode::getClass<PHG4HitContainer>(topNode,
0285                                                                     rec._name);
0286       if (!hits)
0287       {
0288         if (_event < 2)
0289         {
0290           std::cout
0291               << "PHG4DSTReader::process_event - Error - can not find node "
0292               << rec._name << std::endl;
0293         }
0294       }
0295       else
0296       {
0297         PHG4HitContainer::ConstRange hit_range = hits->getHits();
0298
0299         if (Verbosity() >= 2)
0300         {
0301           std::cout << "PHG4DSTReader::process_event - processing "
0302                     << rec._name << " and received " << hits->size() << " hits"
0303                     << std::endl;
0304         }
0305
0306         for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first;
0307              hit_iter != hit_range.second; hit_iter++)
0308         {
0309           PHG4Hit *hit = hit_iter->second;
0310
0311           //                  hit_type * hit = dynamic_cast<hit_type *>(hit_raw);
0312           //
0313           assert(hit);
0314
0315           new ((*(rec._arr.get()))[rec._cnt]) hit_type(hit);
0316
0317           hit_type *new_hit =
0318               dynamic_cast<hit_type *>(rec._arr.get()->At(rec._cnt));
0319           assert(new_hit);
0320
0321           //                  for (int i = 0; i < 2; i++)
0322           //                    {
0323           //                      new_hit->set_x(i, hit->get_x(i));
0324           //                      new_hit->set_y(i, hit->get_y(i));
0325           //                      new_hit->set_z(i, hit->get_z(i));
0326           //                      new_hit->set_t(i, hit->get_t(i));
0327           //                    }
0328           //                  new_hit->set_edep(hit->get_edep());
0329           //                  new_hit->set_layer(hit->get_layer());
0330           //                  new_hit->set_hit_id(hit->get_hit_id());
0331           //                  new_hit->set_trkid(hit->get_trkid());
0332
0333           //                  *new_hit = (*hit);
0334
0335           if (_load_active_particle)
0336           {
0337             _particle_set.insert(hit->get_trkid());
0338           }
0339
0340           if (Verbosity() >= 2)
0341           {
0342             std::cout << "PHG4DSTReader::process_event - processing "
0343                       << rec._name << " and hit " << hit->get_hit_id()
0344                       << " with track id " << hit->get_trkid() << std::endl;
0345           }
0346
0347           rec._cnt++;
0348         }
0349       }  // if (!hits)
0350     }    //      if (rec._type == record::typ_hit)
0351     else if (rec._type == record::typ_tower)
0352     {
0353       if (Verbosity() >= 2)
0354       {
0355         std::cout << "PHG4DSTReader::process_event - processing tower "
0356                   << rec._name << std::endl;
0357       }
0358
0359       RawTowerContainer *hits = findNode::getClass<RawTowerContainer>(
0360           topNode, rec._name);
0361       if (!hits)
0362       {
0363         if (_event < 2)
0364         {
0365           std::cout
0366               << "PHG4DSTReader::process_event - Error - can not find node "
0367               << rec._name << std::endl;
0368         }
0369       }
0370       else
0371       {
0372         RawTowerContainer::ConstRange hit_range = hits->getTowers();
0373
0374         if (Verbosity() >= 2)
0375         {
0376           std::cout << "PHG4DSTReader::process_event - processing "
0377                     << rec._name << " and received " << hits->size()
0378                     << " tower hits" << std::endl;
0379         }
0380
0381         for (RawTowerContainer::ConstIterator hit_iter = hit_range.first;
0382              hit_iter != hit_range.second; hit_iter++)
0383         {
0384           RawTower *hit_raw = hit_iter->second;
0385
0386           RawTower_type *hit = dynamic_cast<RawTower_type *>(hit_raw);
0387           //                  RawTower * hit = hit_iter->second;
0388
0389           assert(hit);
0390
0391           if (hit->get_energy() < _tower_zero_sup)
0392           {
0393             if (Verbosity() >= 2)
0394             {
0395               std::cout
0396                   << "PHG4DSTReader::process_event - suppress low energy tower hit "
0397                   << rec._name << " @ ("
0398                   //                            << hit->get_thetaMin()
0399                   //                            << ", " << hit->get_phiMin()
0400                   << "), Energy = " << hit->get_energy() << std::endl;
0401             }
0402
0403             continue;
0404           }
0405
0406           new ((*(rec._arr.get()))[rec._cnt]) RawTower_type();
0407
0408           if (Verbosity() >= 2)
0409           {
0410             std::cout
0411                 << "PHG4DSTReader::process_event - processing Tower hit "
0412                 << rec._name << " @ ("
0413                 //                        << hit->get_thetaMin() << ", "
0414                 //                        << hit->get_phiMin()
0415                 << "), Energy = " << hit->get_energy() << " - "
0416                 << rec._arr.get()->At(rec._cnt)->ClassName() << std::endl;
0417           }
0418
0419           //                  rec._arr->Print();
0420
0421           RawTower_type *new_hit =
0422               dynamic_cast<RawTower_type *>(rec._arr.get()->At(rec._cnt));
0423           assert(new_hit);
0424
0425           //                  //v1 copy
0426           //                  new_hit->get_bineta(hit->get_bineta());
0427           //                  new_hit->get_binphi(hit->get_binphi());
0428           //
0429           //                  //v2 copy
0430           //                  new_hit->set_thetaMin(hit->get_thetaMin());
0431           //                  new_hit->set_thetaSize(hit->get_thetaSize());
0432           //                  new_hit->set_phiMin(hit->get_phiMin());
0433           //                  new_hit->set_phiSize(hit->get_phiSize());
0434           //                  new_hit->set_zMin(hit->get_zMin());
0435           //                  new_hit->set_zSize(hit->get_zSize());
0436
0437           *new_hit = (*hit);
0438
0439           rec._cnt++;
0440         }
0441       }  // if (!hits)
0442     }    //      if (rec._type == record::typ_hit)
0443     else if (rec._type == record::typ_jets)
0444     {
0445       //          std::std::cout
0446       //              << "PHG4DSTReader::AddJet - Error - temp. disabled until jet added back to sPHENIX software"
0447       //              << std::std::endl;
0448       //
0449       if (Verbosity() >= 2)
0450       {
0451         std::cout << "PHG4DSTReader::process_event - processing jets "
0452                   << rec._name << std::endl;
0453       }
0454
0455       JetContainer *hits = findNode::getClass<JetContainer>(topNode, rec._name);
0456       if (!hits)
0457       {
0458         if (_event < 2)
0459         {
0460           std::cout
0461               << "PHG4DSTReader::process_event - Error - can not find node "
0462               << rec._name << std::endl;
0463         }
0464       }
0465       else
0466       {
0467         if (Verbosity() >= 2)
0468         {
0469           std::cout << "PHG4DSTReader::process_event - processing "
0470                     << rec._name << " and received " << hits->size() << " jets"
0471                     << std::endl;
0472         }
0473
0474         // for every recojet
0475         for (auto hit_raw : *hits)
0476         {
0477           /* Jet *hit_raw = iter.second; */
0478
0479           if (Verbosity() >= 2)
0480           {
0481             std::cout << "PHG4DSTReader::process_event - processing jet "
0482                       << rec._name << " @ (" << hit_raw->get_eta() << ", "
0483                       << hit_raw->get_phi() << "), pT = " << hit_raw->get_pt()
0484                       << " - with raw type " << hit_raw->ClassName() << std::endl;
0485           }
0486
0487           PHPyJet_type *hit = dynamic_cast<PHPyJet_type *>(hit_raw);
0488
0489           assert(hit);
0490
0491           new ((*(rec._arr.get()))[rec._cnt]) PHPyJet_type();
0492
0493           PHPyJet_type *new_hit =
0494               dynamic_cast<PHPyJet_type *>(rec._arr.get()->At(rec._cnt));
0495           assert(new_hit);
0496
0497           *new_hit = (Jetv2) (*hit);
0498
0499           rec._cnt++;
0500         }
0501       }  // if (!hits)
0502     }    //      if (rec._type == record::typ_hit)
0503     else if (rec._type == record::typ_part)
0504     {
0505       if (_load_all_particle)
0506       {
0507         static bool once = true;
0508         if (once)
0509         {
0510           std::cout
0511               << "PHG4DSTReader::process_event - will load all particle from G4TruthInfo"
0512               << std::endl;
0513
0514           once = false;
0515         }
0516
0517         for (auto particle_iter : truthInfoList->GetMap())
0518         {
0519           _particle_set.insert(particle_iter.first);
0520         }
0521
0522       }  //          if (_load_all_particle)
0523       else
0524       {
0525         static bool once = true;
0526         if (once)
0527         {
0528           if (Verbosity() > 0)
0529           {
0530             std::cout << "PHG4DSTReader::process_event - will load primary particle from G4TruthInfo" << std::endl;
0531           }
0532           once = false;
0533         }
0534
0535         PHG4TruthInfoContainer::ConstRange primary_range =
0536             truthInfoList->GetPrimaryParticleRange();
0537
0538         for (PHG4TruthInfoContainer::ConstIterator particle_iter =
0539                  primary_range.first;
0540              particle_iter != primary_range.second;
0541              ++particle_iter)
0542         {
0543           _particle_set.insert(particle_iter->first);
0544
0545         }  //          if (_load_all_particle) else
0546       }
0547       for (int i : _particle_set)
0548       {
0549         auto particle = truthInfoList->GetParticle(i);
0550         if (!particle)
0551         {
0552           std::cout
0553               << "PHG4DSTReader::process_event - ERROR - can not find particle/track ID "
0554               << i << " in G4TruthInfo" << std::endl;
0555
0556           continue;
0557         }
0558
0559         add_particle(rec, particle);
0560       }  // for(PartSet_t::const_iterator i = _particle_set.begin();i!=_particle_set.end();i++)
0561
0562     }  //      if (rec._type == record::typ_part)
0563     else if (rec._type == record::typ_vertex)
0564     {
0565       static bool once = true;
0566       if (once)
0567       {
0568         if (Verbosity() > 0)
0569         {
0570           std::cout << "PHG4DSTReader::process_event - will load vertex from G4TruthInfo" << std::endl;
0571         }
0572         once = false;
0573       }
0574
0575       for (int i : _vertex_set)
0576       {
0577         PHG4VtxPoint *v = truthInfoList->GetVtx(i);
0578         if (!v)
0579         {
0580           std::cout
0581               << "PHG4DSTReader::process_event - ERROR - can not find vertex ID "
0582               << i << " in G4TruthInfo" << std::endl;
0583
0584           continue;
0585         }
0586
0587         new ((*(rec._arr.get()))[rec._cnt]) vertex_type();
0588
0589         if (Verbosity() >= 2)
0590         {
0591           std::cout << "PHG4DSTReader::process_event - saving vertex id "
0592                     << i << std::endl;
0593         }
0594
0595         vertex_type *new_v =
0596             static_cast<vertex_type *>(rec._arr.get()->At(rec._cnt));
0597         assert(new_v);
0598
0599         new_v->set_x(v->get_x());
0600         new_v->set_y(v->get_y());
0601         new_v->set_z(v->get_z());
0602         new_v->set_t(v->get_t());
0603         new_v->set_id(v->get_id());
0604
0605         rec._cnt++;
0606       }  // else if (rec._type == record::typ_vertex)
0607
0608     }  //          if (_load_all_particle)
0609
0610   }  //  for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
0611
0612   if (_T)
0613   {
0614     _T->Fill();
0615   }
0616
0617   return 0;
0618 }  //  for (records_t::iterator it = _records.begin(); it != _records.end(); ++it)
0619
0620 void PHG4DSTReader::add_particle(PHG4DSTReader::record &rec, PHG4Particle *part)
0621 {
0622   assert(part);
0623
0624   //              if (Verbosity() >= 2)
0625   //                std::cout << "PHG4DSTReader::process_event - Particle type is "
0626   //                    << p_raw->GetName() << " - " << p_raw->ClassName()
0627   //                    << " get_track_id = " << p_raw->get_track_id() << std::endl;
0628
0629   //              part_type * part = dynamic_cast<part_type *>(p_raw);
0630
0631   //              assert(part);
0632
0633   new ((*(rec._arr.get()))[rec._cnt]) part_type();
0634
0635   part_type *new_part = static_cast<part_type *>(rec._arr.get()->At(rec._cnt));
0636   assert(new_part);
0637
0638   new_part->set_track_id(part->get_track_id());
0639   new_part->set_vtx_id(part->get_vtx_id());
0640   new_part->set_parent_id(part->get_parent_id());
0641   new_part->set_primary_id(part->get_primary_id());
0642   new_part->set_name(part->get_name());
0643   new_part->set_pid(part->get_pid());
0644   new_part->set_px(part->get_px());
0645   new_part->set_py(part->get_py());
0646   new_part->set_pz(part->get_pz());
0647   new_part->set_e(part->get_e());
0648
0649   _vertex_set.insert(part->get_vtx_id());
0650
0651   rec._cnt++;
0652 }
0653
0654 int PHG4DSTReader::End(PHCompositeNode * /*topNode*/)
0655 {
0656   if (Verbosity() > 1)
0657   {
0658     std::cout << "PHG4DSTReader::End - Clean ups" << std::endl;
0659   }
0660
0661   if (_T)
0662   {
0663     PHTFileServer::get().cd(_out_file_name);
0664     _T->Write();
0665     _T->ResetBranchAddresses();
0666   }
0667
0668   _records.clear();
0669
0670   return 0;
0671 }