QA/Tracking/SiliconSeedsQA.cc

0001
0002 #include "SiliconSeedsQA.h"
0003
0004 #include <fun4all/Fun4AllHistoManager.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006
0007 #include <qautils/QAHistManagerDef.h>
0008 #include <qautils/QAUtil.h>
0009
0010 #include <globalvertex/SvtxVertex.h>
0011 #include <globalvertex/SvtxVertexMap.h>
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/getClass.h>
0014
0015 #include <trackbase/TrackFitUtils.h>
0016 #include <trackbase/TrkrClusterContainer.h>
0017
0018 #include <trackbase_historic/SvtxTrack.h>
0019 #include <trackbase_historic/SvtxTrackMap.h>
0020 #include <trackbase_historic/TrackAnalysisUtils.h>
0021
0022 #include <TH2.h>
0023 #include <TProfile.h>
0024 #include <TProfile2D.h>
0025
0026 #include <trackbase/ActsGeometry.h>
0027 //____________________________________________________________________________..
0028 SiliconSeedsQA::SiliconSeedsQA(const std::string &name)
0029   : SubsysReco(name)
0030 {
0031 }
0032
0033 //____________________________________________________________________________..
0034 int SiliconSeedsQA::InitRun(PHCompositeNode * /*unused*/)
0035 {
0036   createHistos();
0037   return Fun4AllReturnCodes::EVENT_OK;
0038 }
0039
0040 //____________________________________________________________________________..
0041 int SiliconSeedsQA::process_event(PHCompositeNode *topNode)
0042 {
0043   auto *clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
0044   auto *geometry = findNode::getClass<ActsGeometry>(topNode, m_actsgeometryName);
0045   auto *trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0046   auto *vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vertexMapName);
0047
0048   if (!trackmap or !clustermap or !geometry or !vertexmap)
0049   {
0050     std::cout << PHWHERE << "Missing node(s), can't continue" << std::endl;
0051     return Fun4AllReturnCodes::ABORTEVENT;
0052   }
0053
0054   h_ntrack1d->Fill(trackmap->size());
0055
0056   std::pair<int, int> ntrack_isfromvtx;    // first: number of tracks associated to a vertex, second: number of tracks not associated to a vertex
0057   std::pair<int, int> ntrack_isposcharge;  // first: number of tracks with negative charge, second: number of tracks with positive charge
0058
0059   for (const auto &[key, track] : *trackmap)
0060   {
0061     if (!track)
0062     {
0063       continue;
0064     }
0065
0066     auto ckeys = get_cluster_keys(track);
0067     float eta = track->get_eta();
0068     float phi = track->get_phi();
0069     float pt = track->get_pt();
0070     int charge = track->get_charge();
0071     float chi2ndf = track->get_quality();
0072
0073     int trkcrossing = track->get_crossing();
0074     // std::cout << "track crossing: " << trkcrossing << std::endl;
0075
0076     int nmaps = 0;
0077     int nintt = 0;
0078     //    int ntpc = 0;
0079     //    int nmms = 0;
0080
0081     for (auto &ckey : ckeys)
0082     {
0083       switch (TrkrDefs::getTrkrId(ckey))
0084       {
0085       case TrkrDefs::mvtxId:
0086         nmaps++;
0087         break;
0088       case TrkrDefs::inttId:
0089         nintt++;
0090         break;
0091       // case TrkrDefs::tpcId:
0092       //   ntpc++;
0093       //   break;
0094       // case TrkrDefs::micromegasId:
0095       //   nmms++;
0096       //   break;
0097       default:
0098         break;
0099       }
0100     }
0101
0102     Acts::Vector3 zero = Acts::Vector3::Zero();
0103     auto dcapair_origin = TrackAnalysisUtils::get_dca(track, zero);
0104
0105     auto *trackvtx = vertexmap->get(track->get_vertex_id());
0106     if (!trackvtx)
0107     {
0108       ntrack_isfromvtx.first++;
0109       continue;
0110     }
0111     ntrack_isfromvtx.second++;
0112
0113     if (charge > 0)
0114     {
0115       ntrack_isposcharge.second++;
0116     }
0117     else
0118     {
0119       ntrack_isposcharge.first++;
0120     }
0121
0122     Acts::Vector3 track_vtx(trackvtx->get_x(), trackvtx->get_y(), trackvtx->get_z());
0123     auto dcapair_vtx = TrackAnalysisUtils::get_dca(track, track_vtx);
0124
0125     h_ntrack->Fill(eta, phi);
0126     h_nmaps->Fill(nmaps);
0127     h_nintt->Fill(nintt);
0128     h_nmaps_nintt->Fill(nmaps, nintt);
0129     h_avgnclus_eta_phi->Fill(eta, phi, nmaps + nintt);
0130     h_trackcrossing->Fill(trkcrossing);
0131     h_trackchi2ndf->Fill(chi2ndf);
0132     h_dcaxyorigin_phi->Fill(phi, dcapair_origin.first.first);
0133     h_dcaxyvtx_phi->Fill(phi, dcapair_vtx.first.first);
0134     h_dcazorigin_phi->Fill(phi, dcapair_origin.second.first);
0135     h_dcazvtx_phi->Fill(phi, dcapair_vtx.second.first);
0136     h_trackpt_inclusive->Fill(pt);
0137     if (charge > 0)
0138     {
0139       h_trackpt_pos->Fill(pt);
0140     }
0141     else
0142     {
0143       h_trackpt_neg->Fill(pt);
0144     }
0145   }
0146
0147   h_ntrack_isfromvtx->SetBinContent(1, h_ntrack_isfromvtx->GetBinContent(1) + ntrack_isfromvtx.first);
0148   h_ntrack_isfromvtx->SetBinContent(2, h_ntrack_isfromvtx->GetBinContent(2) + ntrack_isfromvtx.second);
0149
0150   h_ntrack_IsPosCharge->SetBinContent(1, h_ntrack_IsPosCharge->GetBinContent(1) + ntrack_isposcharge.first);
0151   h_ntrack_IsPosCharge->SetBinContent(2, h_ntrack_IsPosCharge->GetBinContent(2) + ntrack_isposcharge.second);
0152
0153   // vertex
0154   h_nvertex->Fill(vertexmap->size());
0155   for (const auto &[key, vertex] : *vertexmap)
0156   {
0157     if (!vertex)
0158     {
0159       continue;
0160     }
0161
0162     float vx = vertex->get_x();
0163     float vy = vertex->get_y();
0164     float vz = vertex->get_z();
0165     float vt = vertex->get_t0();
0166     float vchi2 = vertex->get_chisq();
0167     int vndof = vertex->get_ndof();
0168     int vcrossing = vertex->get_beam_crossing();
0169
0170     // std::cout << "vertex (x,y,z,t,chi2,ndof,crossing)=(" << vx << "," << vy << "," << vz << "," << vt << "," << vchi2 << "," << vndof << "," << vcrossing << ")" << std::endl;
0171
0172     h_vx->Fill(vx);
0173     h_vy->Fill(vy);
0174     h_vx_vy->Fill(vx, vy);
0175     h_vz->Fill(vz);
0176     h_vt->Fill(vt);
0177     h_vchi2dof->Fill(float(vchi2 / vndof));
0178     h_vcrossing->Fill(vcrossing);
0179
0180     h_ntrackpervertex->Fill(vertex->size_tracks());
0181   }
0182
0183   return Fun4AllReturnCodes::EVENT_OK;
0184 }
0185 std::vector<TrkrDefs::cluskey> SiliconSeedsQA::get_cluster_keys(SvtxTrack *track)
0186 {
0187   std::vector<TrkrDefs::cluskey> out;
0188   for (const auto &seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0189   {
0190     if (seed)
0191     {
0192       std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0193     }
0194   }
0195   return out;
0196 }
0197
0198 //____________________________________________________________________________..
0199 int SiliconSeedsQA::EndRun(const int /*runnumber*/)
0200 {
0201   return Fun4AllReturnCodes::EVENT_OK;
0202 }
0203
0204 //____________________________________________________________________________..
0205 int SiliconSeedsQA::End(PHCompositeNode * /*unused*/)
0206 {
0207   return Fun4AllReturnCodes::EVENT_OK;
0208 }
0209
0210 std::string SiliconSeedsQA::getHistoPrefix() const
0211 {
0212   return std::string("h_") + Name() + std::string("_");
0213 }
0214
0215 void SiliconSeedsQA::createHistos()
0216 {
0217   auto *hm = QAHistManagerDef::getHistoManager();
0218   assert(hm);
0219
0220   {
0221     h_nmaps = new TH1F(std::string(getHistoPrefix() + "nmaps").c_str(), "MVTX clusters per track;Number of MVTX clusters per track;Entries", 7, -0.5, 6.5);
0222     hm->registerHisto(h_nmaps);
0223   }
0224
0225   {
0226     h_nintt = new TH1F(std::string(getHistoPrefix() + "nintt").c_str(), "INTT clusters per track;Number of INTT clusters per track;Entries", 7, -0.5, 6.5);
0227     hm->registerHisto(h_nintt);
0228   }
0229
0230   {
0231     h_nmaps_nintt = new TH2F(std::string(getHistoPrefix() + "nmaps_nintt").c_str(), "MVTX vs INTT clusters per track;Number of MVTX clusters per track;Number of INTT clusters per track;Entries", 7, -0.5, 6.5, 7, -0.5, 6.5);
0232     hm->registerHisto(h_nmaps_nintt);
0233   }
0234
0235   {
0236     h_ntrack1d = new TH1F(std::string(getHistoPrefix() + "nrecotracks1d").c_str(), "Number of reconstructed tracks;Number of silicon tracklets;Entries", 50, 0, 200);
0237     hm->registerHisto(h_ntrack1d);
0238   }
0239
0240   {
0241     h_ntrack = new TH2F(std::string(getHistoPrefix() + "nrecotracks").c_str(), "Number of reconstructed tracks;#eta;#phi [rad];Entries", 100, -1.1, 1.1, 300, -3.14159, 3.1459);
0242     hm->registerHisto(h_ntrack);
0243   }
0244
0245   {
0246     h_avgnclus_eta_phi = new TProfile2D(std::string(getHistoPrefix() + "avgnclus_eta_phi").c_str(), "Average number of clusters per track;#eta;#phi [rad];Average number of clusters per track", 100, -1.1, 1.1, 300, -3.14159, 3.1459, 0, 10);
0247     hm->registerHisto(h_avgnclus_eta_phi);
0248   }
0249
0250   {
0251     h_trackcrossing = new TH1F(std::string(getHistoPrefix() + "trackcrossing").c_str(), "Track beam bunch crossing;Track crossing;Entries", 110, -10, 100);
0252     hm->registerHisto(h_trackcrossing);
0253   }
0254
0255   {
0256     h_trackchi2ndf = new TH1F(std::string(getHistoPrefix() + "trackchi2ndf").c_str(), "Track chi2/ndof;Track #chi2/ndof;Entries", 100, 0, 20);
0257     hm->registerHisto(h_trackchi2ndf);
0258   }
0259
0260   {
0261     h_dcaxyorigin_phi = new TH2F(std::string(getHistoPrefix() + "dcaxyorigin_phi").c_str(), "DCA xy origin vs phi;#phi [rad];DCA_{xy} wrt origin [cm];Entries", 300, -3.14159, 3.1459, 90, -3, 3);
0262     hm->registerHisto(h_dcaxyorigin_phi);
0263   }
0264
0265   {
0266     h_dcaxyvtx_phi = new TH2F(std::string(getHistoPrefix() + "dcaxyvtx_phi").c_str(), "DCA xy vertex vs phi;#phi [rad];DCA_{xy} wrt vertex [cm];Entries", 300, -3.14159, 3.1459, 90, -3, 3);
0267     hm->registerHisto(h_dcaxyvtx_phi);
0268   }
0269
0270   {
0271     h_dcazorigin_phi = new TH2F(std::string(getHistoPrefix() + "dcazorigin_phi").c_str(), "DCA z origin vs phi;#phi [rad];DCA_{z} wrt origin [cm];Entries", 300, -3.14159, 3.1459, 160, -20, 20);
0272     hm->registerHisto(h_dcazorigin_phi);
0273   }
0274
0275   {
0276     h_dcazvtx_phi = new TH2F(std::string(getHistoPrefix() + "dcazvtx_phi").c_str(), "DCA z vertex vs phi;#phi [rad];DCA_{z} wrt vertex [cm];Entries", 300, -3.14159, 3.1459, 160, -20, 20);
0277     hm->registerHisto(h_dcazvtx_phi);
0278   }
0279
0280   {
0281     h_ntrack_isfromvtx = new TH1F(std::string(getHistoPrefix() + "ntrack_isfromvtx").c_str(), "Num of tracks associated to a vertex;Is track associated to a vertex;Entries", 2, -0.5, 1.5);
0282     hm->registerHisto(h_ntrack_isfromvtx);
0283   }
0284
0285   {
0286     h_trackpt_inclusive = new TH1F(std::string(getHistoPrefix() + "trackpt").c_str(), "Track p_{T};p_{T} [GeV/c];Entries", 100, 0, 10);
0287     hm->registerHisto(h_trackpt_inclusive);
0288   }
0289
0290   {
0291     h_trackpt_pos = new TH1F(std::string(getHistoPrefix() + "trackpt_pos").c_str(), "Track p_{T} positive;Positive track p_{T} [GeV/c];Entries", 100, 0, 10);
0292     hm->registerHisto(h_trackpt_pos);
0293   }
0294
0295   {
0296     h_trackpt_neg = new TH1F(std::string(getHistoPrefix() + "trackpt_neg").c_str(), "Track p_{T} negative;Negative track p_{T} [GeV/c];Entries", 100, 0, 10);
0297     hm->registerHisto(h_trackpt_neg);
0298   }
0299
0300   {
0301     h_ntrack_IsPosCharge = new TH1F(std::string(getHistoPrefix() + "ntrack_IsPosCharge").c_str(), "Num of tracks with positive charge;Is track positive charged;Entries", 2, -0.5, 1.5);
0302     hm->registerHisto(h_ntrack_IsPosCharge);
0303   }
0304
0305   // vertex
0306   {
0307     h_nvertex = new TH1F(std::string(getHistoPrefix() + "nrecovertices").c_str(), "Num of reco vertices per event;Number of vertices;Entries", 20, 0, 20);
0308     hm->registerHisto(h_nvertex);
0309   }
0310
0311   {
0312     h_vx = new TH1F(std::string(getHistoPrefix() + "vx").c_str(), "Vertex x;Vertex x [cm];Entries", 100, -2.5, 2.5);
0313     hm->registerHisto(h_vx);
0314   }
0315
0316   {
0317     h_vy = new TH1F(std::string(getHistoPrefix() + "vy").c_str(), "Vertex y;Vertex y [cm];Entries", 100, -2.5, 2.5);
0318     hm->registerHisto(h_vy);
0319   }
0320
0321   {
0322     h_vx_vy = new TH2F(std::string(getHistoPrefix() + "vx_vy").c_str(), "Vertex x vs y;Vertex x [cm];Vertex y [cm];Entries", 100, -2.5, 2.5, 100, -2.5, 2.5);
0323     hm->registerHisto(h_vx_vy);
0324   }
0325
0326   {
0327     h_vz = new TH1F(std::string(getHistoPrefix() + "vz").c_str(), "Vertex z;Vertex z [cm];Entries", 50, -25, 25);
0328     hm->registerHisto(h_vz);
0329   }
0330
0331   {
0332     h_vt = new TH1F(std::string(getHistoPrefix() + "vt").c_str(), "Vertex t;Vertex t [ns];Entries", 100, -1000, 20000);
0333     hm->registerHisto(h_vt);
0334   }
0335
0336   {
0337     h_vcrossing = new TH1F(std::string(getHistoPrefix() + "vertexcrossing").c_str(), "Vertex beam bunch crossing;Vertex crossing;Entries", 100, -100, 300);
0338     hm->registerHisto(h_vcrossing);
0339   }
0340
0341   {
0342     h_vchi2dof = new TH1F(std::string(getHistoPrefix() + "vertexchi2dof").c_str(), "Vertex chi2/ndof;Vertex #chi2/ndof;Entries", 100, 0, 20);
0343     hm->registerHisto(h_vchi2dof);
0344   }
0345
0346   {
0347     h_ntrackpervertex = new TH1F(std::string(getHistoPrefix() + "ntrackspervertex").c_str(), "Num of tracks per vertex;Number of tracks per vertex;Entries", 50, 0, 50);
0348     hm->registerHisto(h_ntrackpervertex);
0349   }
0350 }