Test/RecoInfoExport/RecoInfoExport.C

0001 #include "RecoInfoExport.h"
0002
0003 #include <phool/getClass.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <phool/PHCompositeNode.h>
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <phool/getClass.h>
0008
0009 #include <g4main/PHG4TruthInfoContainer.h>
0010 #include <g4main/PHG4Particle.h>
0011 #include <g4main/PHG4VtxPoint.h>
0012 #include <g4main/PHG4Hit.h>
0013
0014 #include <g4hough/SvtxVertexMap.h>
0015 #include <g4hough/SvtxVertex.h>
0016 #include <g4hough/SvtxTrackMap.h>
0017 #include <g4hough/SvtxTrack.h>
0018
0019 #include <g4eval/SvtxEvalStack.h>
0020 #include <g4eval/SvtxTrackEval.h>
0021 #include <g4eval/SvtxVertexEval.h>
0022 #include <g4eval/SvtxTruthEval.h>
0023
0024 #include <calobase/RawTowerContainer.h>
0025 #include <calobase/RawTowerGeomContainer.h>
0026 #include <calobase/RawTower.h>
0027 #include <calobase/RawClusterContainer.h>
0028 #include <calobase/RawCluster.h>
0029
0030 #include <TH1D.h>
0031 #include <TH2D.h>
0032 #include <TDatabasePDG.h>
0033 #include <TVector3.h>
0034
0035 #include <iostream>
0036 #include <cassert>
0037 #include <fstream>
0038 #include <sstream>
0039 #include <boost/format.hpp>
0040 #include <boost/math/special_functions/sign.hpp>
0041
0042 using namespace std;
0043
0044 RecoInfoExport::RecoInfoExport(const string &name) :
0045     SubsysReco(name), _event(0), _calo_names(
0046       { "CEMC", "HCALIN", "HCALOUT" }), _tower_threshold(0), _pT_threshold(0), _min_track_hit_dist(
0047         0)
0048 {
0049 }
0050
0051 int
0052 RecoInfoExport::Init(PHCompositeNode *topNode)
0053 {
0054
0055   return 0;
0056 }
0057
0058 int
0059 RecoInfoExport::process_event(PHCompositeNode *topNode)
0060 {
0061   ++_event;
0062
0063   stringstream fname;
0064   fname << _file_prefix << "_Event" << _event << ".dat";
0065   fstream fdata(fname.str(), ios_base::out);
0066
0067   for (auto & calo_name : _calo_names)
0068     {
0069       string towernodename = "TOWER_CALIB_" + calo_name;
0070       // Grab the towers
0071       RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode,
0072           towernodename.c_str());
0073       if (!towers)
0074         {
0075           std::cout << PHWHERE << ": Could not find node "
0076               << towernodename.c_str() << std::endl;
0077           return Fun4AllReturnCodes::ABORTRUN;
0078         }
0079       string towergeomnodename = "TOWERGEOM_" + calo_name;
0080       RawTowerGeomContainer *towergeom = findNode::getClass<
0081           RawTowerGeomContainer>(topNode, towergeomnodename.c_str());
0082       if (!towergeom)
0083         {
0084           cout << PHWHERE << ": Could not find node "
0085               << towergeomnodename.c_str() << endl;
0086           return Fun4AllReturnCodes::ABORTRUN;
0087         }
0088
0089       set<const RawTower *> good_towers;
0090
0091       RawTowerContainer::ConstRange begin_end = towers->getTowers();
0092       RawTowerContainer::ConstIterator rtiter;
0093       for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0094         {
0095           const RawTower *tower = rtiter->second;
0096           assert(tower);
0097           if (tower->get_energy() > _tower_threshold)
0098             {
0099               good_towers.insert(tower);
0100             }
0101         }
0102
0103       fdata
0104           << (boost::format("%1% (1..%2% hits)") % calo_name
0105               % good_towers.size()) << endl;
0106
0107       bool first = true;
0108       for (const auto & tower : good_towers)
0109         {
0110           assert(tower);
0111
0112           float eta = towergeom->get_etacenter(tower->get_bineta());
0113           float phi = towergeom->get_phicenter(tower->get_binphi());
0114
0115           phi = atan2(cos(phi), sin(phi));
0116
0117           if (first)
0118             {
0119               first = false;
0120             }
0121           else
0122             fdata << ",";
0123
0124           fdata
0125               << (boost::format("[%1%,%2%,%3%]") % eta % phi
0126                   % tower->get_energy());
0127
0128         }
0129
0130       fdata << endl;
0131     }
0132
0133     {
0134       fdata << "Track list" << endl;
0135
0136       // need things off of the DST...
0137       PHG4TruthInfoContainer* truthinfo = findNode::getClass<
0138           PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0139       if (!truthinfo)
0140         {
0141           cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
0142           exit(-1);
0143         }
0144
0145       // create SVTX eval stack
0146       SvtxEvalStack svtxevalstack(topNode);
0147
0148 //  SvtxVertexEval* vertexeval = svtxevalstack.get_vertex_eval();
0149 //  SvtxTrackEval* trackeval = svtxevalstack.get_track_eval();
0150       SvtxTruthEval* trutheval = svtxevalstack.get_truth_eval();
0151
0152       // loop over all truth particles
0153       PHG4TruthInfoContainer::Range range =
0154           truthinfo->GetPrimaryParticleRange();
0155       for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0156           iter != range.second; ++iter)
0157         {
0158           PHG4Particle* g4particle = iter->second;
0159
0160           const TVector3 mom(g4particle->get_px(), g4particle->get_py(),
0161               g4particle->get_pz());
0162
0163           std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
0164
0165           map<float, PHG4Hit *> time_sort;
0166           map<float, PHG4Hit *> layer_sort;
0167           for (auto & hit : g4hits)
0168             {
0169               if (hit)
0170                 {
0171                   time_sort[hit->get_avg_t()] = hit;
0172                 }
0173             }
0174
0175           for (auto & hit_pair : time_sort)
0176             {
0177
0178               if (hit_pair.second->get_layer() != UINT_MAX
0179                   and layer_sort.find(hit_pair.second->get_layer())
0180                       == layer_sort.end())
0181                 {
0182                   layer_sort[hit_pair.second->get_layer()] = hit_pair.second;
0183                 }
0184             }
0185
0186           if (layer_sort.size() > 5 and mom.Pt() > _pT_threshold) // minimal track length cut
0187             {
0188
0189               stringstream spts;
0190
0191               TVector3 last_pos(0, 0, 0);
0192
0193               bool first = true;
0194               for (auto & hit_pair : layer_sort)
0195                 {
0196                   TVector3 pos(hit_pair.second->get_avg_x(),
0197                       hit_pair.second->get_avg_y(),
0198                       hit_pair.second->get_avg_z());
0199
0200                   // hit step cuts
0201                   if ((pos - last_pos).Mag() < _min_track_hit_dist
0202                       and hit_pair.first != (layer_sort.rbegin()->first)
0203                       and hit_pair.first != (layer_sort.begin()->first))
0204                     continue;
0205
0206                   last_pos = pos;
0207
0208                   if (first)
0209                     {
0210                       first = false;
0211                     }
0212                   else
0213                     spts << ",";
0214
0215                   spts << "[";
0216                   spts << pos.x();
0217                   spts << ",";
0218                   spts << pos.y();
0219                   spts << ",";
0220                   spts << pos.z();
0221                   spts << "]";
0222                 }
0223
0224               const int abs_pid = abs(g4particle->get_pid());
0225               int t = 5;
0226               if (abs_pid
0227                   == TDatabasePDG::Instance()->GetParticle("pi+")->PdgCode())
0228                 {
0229                   t = 1;
0230                 }
0231               else if (abs_pid
0232                   == TDatabasePDG::Instance()->GetParticle("proton")->PdgCode())
0233                 {
0234                   t = 2;
0235                 }
0236               else if (abs_pid
0237                   == TDatabasePDG::Instance()->GetParticle("K+")->PdgCode())
0238                 {
0239                   t = 3;
0240                 }
0241               else if (abs_pid
0242                   == TDatabasePDG::Instance()->GetParticle("e-")->PdgCode())
0243                 {
0244                   t = 3;
0245                 }
0246
0247               const TParticlePDG * pdg_part =
0248                   TDatabasePDG::Instance()->GetParticle(11);
0249               const int c =
0250                   (pdg_part != nullptr) ? (copysign(1, pdg_part->Charge())) : 0;
0251
0252               fdata
0253                   << (boost::format(
0254                       "{ \"pt\": %1%, \"t\": %2%, \"e\": %3%, \"p\": %4%, \"c\": %5%, \"pts\":[ %6% ]}")
0255                       % mom.Pt() % t % mom.PseudoRapidity() % mom.Phi() % c
0256                       % spts.str()) << endl;
0257
0258             }
0259         }
0260     }
0261
0262   fdata.close();
0263   return 0;
0264 }
0265
0266 int
0267 RecoInfoExport::End(PHCompositeNode *topNode)
0268 {
0269   return 0;
0270 }