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