Tracking/dEdx/dEdxFitter.cc

0001 #include "dEdxFitter.h"
0002
0003 #include <phool/phool.h>
0004 #include <phool/getClass.h>
0005 #include <fun4all/PHTFileServer.h>
0006 #include <fun4all/Fun4AllServer.h>
0007 #include <trackbase_historic/TrackAnalysisUtils.h>
0008 #include <g4detectors/PHG4TpcGeom.h>
0009
0010 #include <iostream>
0011
0012 //____________________________________
0013 dEdxFitter::dEdxFitter(const std::string &name):
0014     SubsysReco(name)
0015 {
0016   //initialize
0017   fitter = std::make_unique<GlobaldEdxFitter>();
0018 }
0019
0020 //___________________________________
0021 int dEdxFitter::InitRun(PHCompositeNode *topNode)
0022 {
0023   std::cout << PHWHERE << " Opening file " << _outfile << std::endl;
0024   PHTFileServer::get().open( _outfile, "RECREATE");
0025
0026   return 0;
0027 }
0028
0029 //__________________________________
0030 //Call user instructions for every event
0031 int dEdxFitter::process_event(PHCompositeNode *topNode)
0032 {
0033   _event++;
0034   if(_event%1000==0) std::cout << PHWHERE << "Events processed: " << _event << std::endl;
0035
0036   GetNodes(topNode);
0037
0038   if(Verbosity()>1)
0039   {
0040     std::cout << "--------------------------------" << std::endl;
0041     std::cout << "event " << _event << std::endl;
0042   }
0043
0044   process_tracks(topNode);
0045
0046   return 0;
0047 }
0048
0049 //_____________________________________
0050 void dEdxFitter::process_tracks(PHCompositeNode *topNode)
0051 {
0052
0053   for(const auto &[key, track] : *_trackmap)
0054   {
0055     if(!track) continue;
0056
0057     double trackID = track->get_id();
0058     if(Verbosity()>1) std::cout << "track ID " << trackID << std::endl;
0059     if(std::isnan(track->get_x()) ||
0060        std::isnan(track->get_y()) ||
0061        std::isnan(track->get_z()) ||
0062        std::isnan(track->get_px()) ||
0063        std::isnan(track->get_py()) ||
0064        std::isnan(track->get_pz()))
0065     {
0066       std::cout << "malformed track:" << std::endl;
0067       track->identify();
0068       std::cout << "skipping..." << std::endl;
0069       continue;
0070     }
0071
0072     // ignore TPC-only tracks
0073     if(!track->get_silicon_seed())
0074     {
0075       if(Verbosity()>1) std::cout << "TPC-only track, skipping..." << std::endl;
0076       continue;
0077     }
0078
0079     std::tuple<int,int,int> nclus = get_nclus(track);
0080     int nmaps = std::get<0>(nclus);
0081     int nintt = std::get<1>(nclus);
0082     int ntpc = std::get<2>(nclus);
0083
0084     if(nmaps>=nmaps_cut && nintt>=nintt_cut && ntpc>=ntpc_cut && fabs(track->get_eta())<eta_cut && get_dcaxy(track)<dcaxy_cut)
0085     {
0086       fitter->addTrack(get_dedx(track),track->get_p());
0087     }
0088
0089     if(fitter->getNtracks() > ntracks_to_fit)
0090     {
0091       minima.push_back(fitter->get_minimum());
0092       fitter->reset();
0093     }
0094   }
0095 }
0096
0097 std::tuple<int,int,int> dEdxFitter::get_nclus(SvtxTrack* track)
0098 {
0099   int nmaps = 0;
0100   int nintt = 0;
0101   int ntpc = 0;
0102
0103   for(auto it = track->get_silicon_seed()->begin_cluster_keys(); it != track->get_silicon_seed()->end_cluster_keys(); ++it)
0104   {
0105     TrkrDefs::cluskey ckey = *it;
0106     auto trkrid = TrkrDefs::getTrkrId(ckey);
0107     if(trkrid == TrkrDefs::mvtxId)
0108     {
0109       nmaps++;
0110     }
0111     else if(trkrid == TrkrDefs::inttId)
0112     {
0113       nintt++;
0114     }
0115   }
0116   for(auto it = track->get_tpc_seed()->begin_cluster_keys(); it != track->get_tpc_seed()->end_cluster_keys(); ++it)
0117   {
0118     ntpc++;
0119   }
0120
0121   return std::make_tuple(nmaps,nintt,ntpc);
0122 }
0123
0124 double dEdxFitter::get_dedx(SvtxTrack* track)
0125 {
0126   float layerThicknesses[4] = {0.0, 0.0, 0.0, 0.0};
0127   // These are randomly chosen layer thicknesses for the TPC, to get the
0128   // correct region thicknesses in an easy to pass way to the helper fxn
0129   layerThicknesses[0] = _tpcgeom->GetLayerCellGeom(7)->get_thickness();
0130   layerThicknesses[1] = _tpcgeom->GetLayerCellGeom(8)->get_thickness();
0131   layerThicknesses[2] = _tpcgeom->GetLayerCellGeom(27)->get_thickness();
0132   layerThicknesses[3] = _tpcgeom->GetLayerCellGeom(50)->get_thickness();
0133
0134   return TrackAnalysisUtils::calc_dedx(track->get_tpc_seed(), _clustermap, _geometry, layerThicknesses);
0135 }
0136
0137 double dEdxFitter::get_dcaxy(SvtxTrack* track)
0138 {
0139   auto vertexit = _vertexmap->find(track->get_vertex_id());
0140   if(vertexit != _vertexmap->end())
0141   {
0142     SvtxVertex* vtx = vertexit->second;
0143     Acts::Vector3 vertex(vtx->get_x(),vtx->get_y(),vtx->get_z());
0144     auto dcapair = TrackAnalysisUtils::get_dca(track,vertex);
0145     return dcapair.first.first;
0146   }
0147   else
0148   {
0149     return std::numeric_limits<float>::quiet_NaN();
0150   }
0151 }
0152
0153 //___________________________________
0154 void dEdxFitter::GetNodes(PHCompositeNode *topNode)
0155 {
0156
0157   _trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0158   if(!_trackmap && _event<2)
0159   {
0160     std::cout << PHWHERE << " cannot find SvtxTrackMap" << std::endl;
0161   }
0162
0163   _clustermap = findNode::getClass<TrkrClusterContainer>(topNode,"TRKR_CLUSTER");
0164   if(!_clustermap && _event<2)
0165   {
0166     std::cout << PHWHERE << " cannot find TrkrClusterContainer TRKR_CLUSTER" << std::endl;
0167   }
0168
0169   _geometry = findNode::getClass<ActsGeometry>(topNode,"ActsGeometry");
0170   if(!_geometry && _event<2)
0171   {
0172     std::cout << PHWHERE << " cannot find ActsGeometry" << std::endl;
0173   }
0174
0175   _tpcgeom = findNode::getClass<PHG4TpcGeomContainer>(topNode,"TPCGEOMCONTAINER");
0176   if(!_tpcgeom && _event<2)
0177   {
0178     std::cout << PHWHERE << " cannot find PHG4TpcGeomContainer TPCGEOMCONTAINER" << std::endl;
0179   }
0180
0181   _vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0182   if(!_vertexmap && _event<2)
0183   {
0184     std::cout << PHWHERE << " cannot find SvtxVertexMap" << std::endl;
0185   }
0186 }
0187
0188 //______________________________________
0189 int dEdxFitter::End(PHCompositeNode *topNode)
0190 {
0191   if(minima.size()==0)
0192   {
0193     minima.push_back(fitter->get_minimum());
0194   }
0195
0196   PHTFileServer::get().cd( _outfile );
0197
0198   double avg_minimum = 0.;
0199   for(double m : minima)
0200   {
0201     avg_minimum += m;
0202   }
0203   avg_minimum /= (double)minima.size();
0204
0205   TF1* pi_band = new TF1("pi_band","bethe_bloch_new_1D(fabs(x)/[1],[0])");
0206   pi_band->SetParameter(0,avg_minimum);
0207   pi_band->SetParameter(1,dedx_constants::m_pi);
0208   pi_band->Write();
0209
0210   TF1* K_band = new TF1("K_band","bethe_bloch_new_1D(fabs(x)/[1],[0])");
0211   K_band->SetParameter(0,avg_minimum);
0212   K_band->SetParameter(1,dedx_constants::m_K);
0213   K_band->Write();
0214
0215   TF1* p_band = new TF1("p_band","bethe_bloch_new_1D(fabs(x)/[1],[0])");
0216   p_band->SetParameter(0,avg_minimum);
0217   p_band->SetParameter(1,dedx_constants::m_p);
0218   p_band->Write();
0219
0220   TF1* d_band = new TF1("d_band","bethe_bloch_new_1D(fabs(x)/[1],[0])");
0221   d_band->SetParameter(0,avg_minimum);
0222   d_band->SetParameter(1,dedx_constants::m_d);
0223   d_band->Write();
0224
0225   if(Verbosity()>0) std::cout << "dEdxFitter extracted minimum: " << avg_minimum << std::endl;
0226
0227   return 0;
0228 }