File indexing completed on 2025-08-05 08:15:41
0001
0002
0003
0004
0005
0006 #include "TrackToCalo.h"
0007
0008 #include <calobase/RawClusterContainer.h>
0009 #include <calobase/RawTowerGeomContainer.h>
0010 #include <calobase/RawCluster.h>
0011 #include <calobase/RawClusterUtility.h>
0012 #include <calobase/RawTowerDefs.h>
0013 #include <calobase/RawTowerGeom.h>
0014 #include <calobase/TowerInfoContainer.h>
0015 #include <calobase/TowerInfo.h>
0016 #include <calobase/TowerInfoDefs.h>
0017
0018 #include <fun4all/Fun4AllReturnCodes.h>
0019
0020 #include <globalvertex/GlobalVertex.h>
0021 #include <globalvertex/GlobalVertexMap.h>
0022 #include <globalvertex/MbdVertexMap.h>
0023 #include <globalvertex/MbdVertex.h>
0024 #include <globalvertex/SvtxVertexMap.h>
0025 #include <globalvertex/SvtxVertex.h>
0026
0027 #include <phool/getClass.h>
0028 #include <phool/PHCompositeNode.h>
0029
0030 #include <trackbase/TrkrCluster.h>
0031 #include <trackbase/TrkrClusterContainer.h>
0032 #include <trackbase/TrkrClusterCrossingAssocv1.h>
0033 #include <trackbase/TrkrDefs.h>
0034 #include <trackbase/TrkrHitSetContainer.h>
0035 #include <trackbase/TrkrHitSet.h>
0036 #include <trackbase_historic/SvtxTrackMap.h>
0037 #include <trackbase_historic/SvtxTrackState_v1.h>
0038 #include <trackbase_historic/TrackSeedContainer.h>
0039 #include <trackbase_historic/TrackSeed.h>
0040 #include <trackbase_historic/TrackAnalysisUtils.h>
0041 #include <trackreco/ActsPropagator.h>
0042
0043 #include <Acts/Geometry/GeometryIdentifier.hpp>
0044 #include <Acts/MagneticField/ConstantBField.hpp>
0045 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0046 #include <Acts/Surfaces/CylinderSurface.hpp>
0047 #include <Acts/Surfaces/PerigeeSurface.hpp>
0048 #include <Acts/Geometry/TrackingGeometry.hpp>
0049
0050 #include <CLHEP/Vector/ThreeVector.h>
0051 #include <cmath>
0052 #include <vector>
0053
0054 #include <TFile.h>
0055 #include <TTree.h>
0056 #include <TH1F.h>
0057 #include <TH2F.h>
0058
0059
0060 TrackToCalo::TrackToCalo(const std::string &name, const std::string &file):
0061 SubsysReco(name),
0062 _outfilename(file),
0063 _outfile(nullptr),
0064 _tree(nullptr)
0065 {
0066 std::cout << "TrackToCalo::TrackToCalo(const std::string &name, const std::string &file) Calling ctor" << std::endl;
0067 }
0068
0069
0070 TrackToCalo::~TrackToCalo()
0071 {
0072 std::cout << "TrackToCalo::~TrackToCalo() Calling dtor" << std::endl;
0073 }
0074
0075
0076 int TrackToCalo::Init(PHCompositeNode *topNode)
0077 {
0078 std::cout << topNode << std::endl;
0079 std::cout << "TrackToCalo::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0080 delete _outfile;
0081 _outfile = new TFile(_outfilename.c_str(), "RECREATE");
0082 delete _tree;
0083 _tree = new TTree("tree", "A tree with track/calo info");
0084 _tree->Branch("calo_z",&m_calo_z,"calo_z/F");
0085 _tree->Branch("calo_phi",&m_calo_phi,"calo_phi/F");
0086 _tree->Branch("track_z",&m_track_z,"track_z/F");
0087 _tree->Branch("track_phi",&m_track_phi,"track_phi/F");
0088 _tree->Branch("dphi",&m_dphi,"dphi/F");
0089 _tree->Branch("dz",&m_dz,"dz/F");
0090 cnt=0;
0091 return Fun4AllReturnCodes::EVENT_OK;
0092 }
0093
0094
0095 int TrackToCalo::process_event(PHCompositeNode *topNode)
0096 {
0097 std::cout<<"TrackToCalo::process_event event "<<cnt<<std::endl;
0098 cnt++;
0099 ResetTreeVectors();
0100
0101 if(!svtxTrackMap)
0102 {
0103 svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0104 if(!svtxTrackMap)
0105 {
0106 std::cout << "TrackToCalo::process_event : SvtxTrackMap (SvtxTrackMap) not found! Aborting!" << std::endl;
0107 return Fun4AllReturnCodes::ABORTEVENT;
0108 }
0109 }
0110
0111 if ( !rawClusterContainer )
0112 {
0113 rawClusterContainer = findNode::getClass<RawClusterContainer>(topNode, m_RawClusCont_EM_name);
0114 if (!rawClusterContainer)
0115 {
0116 std::cout << "TrackToCalo::process_event : RawClusterContainer (" << m_RawClusCont_EM_name << ") not found! Aborting!" << std::endl;
0117 return Fun4AllReturnCodes::ABORTEVENT;
0118 }
0119 }
0120
0121 if(!trkrClusterContainer)
0122 {
0123 trkrClusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0124 if(!trkrClusterContainer)
0125 {
0126 std::cout << "TrackToCalo::process_event : TrkrClusterContainer (TRKR_CLUSTER) not found! Aborting!" << std::endl;
0127 return Fun4AllReturnCodes::ABORTEVENT;
0128 }
0129 }
0130
0131 if(!rawTowerGeomContainer)
0132 {
0133 rawTowerGeomContainer = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0134 if(!rawTowerGeomContainer)
0135 {
0136 std::cout << "TrackToCalo::process_event : RawTowerGeomContainer (TOWERGEOM_CEMC) not found! Aborting!" << std::endl;
0137 return Fun4AllReturnCodes::ABORTEVENT;
0138 }
0139 }
0140
0141 double caloRadiusEMCal;
0142 if (m_use_emcal_radius)
0143 {
0144 caloRadiusEMCal = m_emcal_radius_user;
0145 }
0146 else
0147 {
0148 caloRadiusEMCal = rawTowerGeomContainer->get_radius();
0149 }
0150
0151 for (auto &iter : *svtxTrackMap)
0152 {
0153 track = iter.second;
0154
0155 if(!track)
0156 {
0157 continue;
0158 }
0159
0160 if(track->get_pt() < m_track_pt_low_cut)
0161 {
0162 continue;
0163 }
0164
0165 int n_tpc_clusters = 0;
0166 tpc_seed = track->get_tpc_seed();
0167 if(tpc_seed)
0168 {
0169 for(auto key_iter = tpc_seed->begin_cluster_keys(); key_iter != tpc_seed->end_cluster_keys(); ++key_iter)
0170 {
0171 const auto& cluster_key = *key_iter;
0172 trkrCluster = trkrClusterContainer->findCluster(cluster_key);
0173 if(!trkrCluster)
0174 {
0175 continue;
0176 }
0177 if(TrkrDefs::getTrkrId(cluster_key) == TrkrDefs::TrkrId::tpcId)
0178 {
0179 n_tpc_clusters++;
0180 }
0181 }
0182 }
0183
0184 if(n_tpc_clusters<m_ntpc_low_cut)
0185 {
0186 continue;
0187 }
0188
0189
0190 thisState = track->get_state(caloRadiusEMCal);
0191
0192 if(!thisState)
0193 {
0194 _track_phi_emc.push_back(NAN);
0195 _track_z_emc.push_back(NAN);
0196 }
0197 else
0198 {
0199 _track_phi_emc.push_back(std::atan2(thisState->get_y(), thisState->get_x()));
0200 _track_z_emc.push_back(thisState->get_z());
0201 }
0202 }
0203
0204 RawClusterContainer::Range begin_end_EMC = rawClusterContainer->getClusters();
0205 RawClusterContainer::Iterator clusIter_EMC;
0206
0207
0208 for (clusIter_EMC = begin_end_EMC.first; clusIter_EMC != begin_end_EMC.second; ++clusIter_EMC)
0209 {
0210 cluster = clusIter_EMC->second;
0211
0212 if(cluster->get_energy() < m_emcal_e_low_cut)
0213 {
0214 continue;
0215 }
0216
0217 _emcal_phi.push_back( std::atan2(cluster->get_y(), cluster->get_x()) );
0218 radius_scale = caloRadiusEMCal / sqrt( pow(cluster->get_x(),2) + pow(cluster->get_y(),2) );
0219 _emcal_z.push_back( radius_scale*cluster->get_z() );
0220 }
0221
0222 for(unsigned int itrack = 0; itrack < _track_phi_emc.size(); itrack++)
0223 {
0224 if(std::isnan(_track_phi_emc.at(itrack)) || std::isnan(_track_z_emc.at(itrack)))
0225 {
0226 continue;
0227 }
0228
0229 for(unsigned int iem = 0; iem < _emcal_z.size(); iem++)
0230 {
0231 m_calo_z = _emcal_z.at(iem);
0232 m_track_z = _track_z_emc.at(itrack);
0233 m_calo_phi = _emcal_phi.at(iem);
0234 m_track_phi = _track_phi_emc.at(itrack);
0235 m_dphi = PiRange(_track_phi_emc.at(itrack) - _emcal_phi.at(iem));
0236 m_dz = _track_z_emc.at(itrack) - _emcal_z.at(iem);
0237 _tree->Fill();
0238 }
0239 }
0240
0241 return Fun4AllReturnCodes::EVENT_OK;
0242 }
0243
0244
0245 int TrackToCalo::End(PHCompositeNode *topNode)
0246 {
0247 std::cout << topNode << std::endl;
0248 _outfile->cd();
0249 _outfile->Write();
0250 _outfile->Close();
0251 return Fun4AllReturnCodes::EVENT_OK;
0252 }
0253
0254 void TrackToCalo::ResetTreeVectors()
0255 {
0256 _track_phi_emc.clear();
0257 _track_z_emc.clear();
0258 _emcal_phi.clear();
0259 _emcal_z.clear();
0260 }
0261
0262 float TrackToCalo::PiRange(float deltaPhi)
0263 {
0264 if(deltaPhi > M_PI)
0265 {
0266 deltaPhi -= 2*M_PI;
0267 }
0268 if(deltaPhi < -M_PI)
0269 {
0270 deltaPhi += 2*M_PI;
0271 }
0272
0273 return deltaPhi;
0274 }
