File indexing completed on 2025-08-06 08:18:50
0001 #include "TpcSeedsQA.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004 #include <qautils/QAUtil.h>
0005
0006 #include <globalvertex/SvtxVertex.h>
0007 #include <globalvertex/SvtxVertexMap.h>
0008
0009 #include <trackbase/ActsGeometry.h>
0010 #include <trackbase/TrackFitUtils.h>
0011 #include <trackbase/TrkrCluster.h>
0012 #include <trackbase/TrkrClusterContainer.h>
0013
0014 #include <g4detectors/PHG4TpcCylinderGeom.h>
0015 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0016
0017 #include <trackbase_historic/SvtxTrack.h>
0018 #include <trackbase_historic/SvtxTrackMap.h>
0019 #include <trackbase_historic/TrackAnalysisUtils.h>
0020 #include <trackbase_historic/TrackSeed.h>
0021 #include <trackbase_historic/TrackSeedContainer.h>
0022
0023 #include <tpc/TpcDistortionCorrectionContainer.h>
0024 #include <tpc/TpcGlobalPositionWrapper.h>
0025
0026 #include <ffarawobjects/Gl1Packet.h>
0027 #include <ffarawobjects/Gl1RawHit.h>
0028 #include <fun4all/Fun4AllHistoManager.h>
0029 #include <fun4all/Fun4AllReturnCodes.h>
0030
0031 #include <phool/PHCompositeNode.h>
0032 #include <phool/getClass.h>
0033
0034 #include <TH2.h>
0035 #include <TH2F.h>
0036 #include <TNtuple.h>
0037 #include <TProfile.h>
0038 #include <TProfile2D.h>
0039
0040 #include <algorithm>
0041 #include <cmath>
0042 #include <format>
0043
0044
0045 TpcSeedsQA::TpcSeedsQA(const std::string &name)
0046 : SubsysReco(name)
0047 {
0048 }
0049
0050
0051 int TpcSeedsQA::InitRun(PHCompositeNode *topNode)
0052 {
0053 createHistos();
0054
0055 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
0056 actsgeom = findNode::getClass<ActsGeometry>(topNode, m_actsGeomName);
0057 g4geom = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, m_g4GeomName);
0058 trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0059 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vertexMapName);
0060
0061 if (!trackmap or !clustermap or !actsgeom or !vertexmap)
0062 {
0063 std::cout << PHWHERE << "Missing node(s), can't continue" << std::endl;
0064 return Fun4AllReturnCodes::ABORTEVENT;
0065 }
0066
0067 if (!g4geom)
0068 {
0069 std::cout << PHWHERE << " unable to find DST node CYLINDERCELLGEOM_SVTX" << std::endl;
0070 return Fun4AllReturnCodes::ABORTRUN;
0071 }
0072
0073
0074 m_globalPositionWrapper.loadNodes(topNode);
0075
0076 m_clusterMover.initialize_geometry(g4geom);
0077 m_clusterMover.set_verbosity(0);
0078
0079 auto *hm = QAHistManagerDef::getHistoManager();
0080 assert(hm);
0081
0082
0083 h_ntrack1d = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks1d")));
0084 h_ntrack1d_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks1d_pos")));
0085 h_ntrack1d_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks1d_neg")));
0086 h_ntrack1d_ptg1 = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks1d_ptg1")));
0087 h_ntrack1d_ptg1_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks1d_ptg1_pos")));
0088 h_ntrack1d_ptg1_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks1d_ptg1_neg")));
0089 h_pt = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "pt")));
0090 h_pt_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "pt_pos")));
0091 h_pt_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "pt_neg")));
0092 h_ntrack_pos = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks_pos")));
0093 h_ntrack_neg = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecotracks_neg")));
0094
0095 h_ntpc_fullpt_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpc_fullpt_pos")));
0096 h_ntpc_fullpt_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpc_fullpt_neg")));
0097 h_ntpc_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpc_pos")));
0098 h_ntpc_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpc_neg")));
0099 h_ntpc_quality_pos = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpc_quality_pos")));
0100 h_ntpc_quality_neg = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpc_quality_neg")));
0101 h_ntpot_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpot_pos")));
0102 h_ntpot_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntpot_neg")));
0103 h_avgnclus_eta_phi_pos = dynamic_cast<TProfile2D *>(hm->getHisto(std::string(getHistoPrefix() + "avgnclus_eta_phi_pos")));
0104 h_avgnclus_eta_phi_neg = dynamic_cast<TProfile2D *>(hm->getHisto(std::string(getHistoPrefix() + "avgnclus_eta_phi_neg")));
0105
0106
0107 h_dcaxyorigin_phi_north_pos = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcaxyorigin_phi_north_pos")));
0108 h_dcaxyorigin_phi_south_pos = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcaxyorigin_phi_south_pos")));
0109 h_dcaxyorigin_phi_north_neg = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcaxyorigin_phi_north_neg")));
0110 h_dcaxyorigin_phi_south_neg = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcaxyorigin_phi_south_neg")));
0111 h_dcaxyvtx_phi_pos = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcaxyvtx_phi_pos")));
0112 h_dcaxyvtx_phi_neg = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcaxyvtx_phi_neg")));
0113 h_dcazorigin_phi_pos = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcazorigin_phi_pos")));
0114 h_dcazorigin_phi_neg = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcazorigin_phi_neg")));
0115 h_dcazvtx_phi_pos = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcazvtx_phi_pos")));
0116 h_dcazvtx_phi_neg = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dcazvtx_phi_neg")));
0117 h_ntrack_isfromvtx_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntrack_isfromvtx_pos")));
0118 h_ntrack_isfromvtx_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntrack_isfromvtx_neg")));
0119 h_cluster_phisize1_fraction_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "cluster_phisize1_fraction_pos")));
0120 h_cluster_phisize1_fraction_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "cluster_phisize1_fraction_neg")));
0121
0122
0123 h_nvertex = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nrecovertices")));
0124 h_vx = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vx")));
0125 h_vy = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vy")));
0126 h_vx_vy = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "vx_vy")));
0127 h_vz = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vz")));
0128 h_vt = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vt")));
0129
0130 h_vchi2dof = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vertexchi2dof")));
0131 h_ntrackpervertex = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "ntrackspervertex")));
0132 h_dedx = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dedx")));
0133 h_mip_dedx = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "mip_dedx")));
0134 h_dedx_pcaz = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "dedx_pcaz")));
0135
0136 nt_sector_event_summary = dynamic_cast<TNtuple *>(hm->getHisto(std::string(getHistoPrefix() + "sector_event_summary")));
0137
0138
0139 std::vector<int> region_layer_low = {7, 23, 39};
0140 std::vector<int> region_layer_high = {22, 38, 54};
0141
0142
0143 const auto range = g4geom->get_begin_end();
0144 for (auto iter = range.first; iter != range.second; ++iter)
0145 {
0146 m_layers.insert(iter->first);
0147
0148 for (int region = 0; region < 3; ++region)
0149 {
0150 if (iter->first >= region_layer_low[region] && iter->first <= region_layer_high[region])
0151 {
0152 m_layerRegionMap.insert(std::make_pair(iter->first, region));
0153 }
0154 }
0155 }
0156
0157 for (const auto ®ion : {0, 1, 2})
0158 {
0159 PhiHistoList phihist;
0160
0161 phihist.cphisize1pT_side0 = h_clusphisize1pt_side0[region];
0162 phihist.cphisize1pT_side1 = h_clusphisize1pt_side1[region];
0163
0164 phihist.cphisizegeq1pT_side0 = h_clusphisizegeq1pt_side0[region];
0165 phihist.cphisizegeq1pT_side1 = h_clusphisizegeq1pt_side1[region];
0166
0167 phihist.Clear();
0168
0169 phihistos.insert(std::make_pair(region, phihist));
0170 }
0171
0172 return Fun4AllReturnCodes::EVENT_OK;
0173 }
0174
0175 float TpcSeedsQA::calc_dedx(TrackSeed *tpcseed)
0176 {
0177 std::vector<TrkrDefs::cluskey> clusterKeys;
0178 clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(),
0179 tpcseed->end_cluster_keys());
0180
0181 std::vector<float> dedxlist;
0182 for (unsigned long cluster_key : clusterKeys)
0183 {
0184 auto detid = TrkrDefs::getTrkrId(cluster_key);
0185 if (detid != TrkrDefs::TrkrId::tpcId)
0186 {
0187 continue;
0188 }
0189 unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
0190 TrkrCluster *cluster = clustermap->findCluster(cluster_key);
0191 float adc = cluster->getAdc();
0192 PHG4TpcCylinderGeom *GeoLayer_local = g4geom->GetLayerCellGeom(layer_local);
0193 float thick = GeoLayer_local->get_thickness();
0194 float r = GeoLayer_local->get_radius();
0195 float alpha = (r * r) / (2 * r * TMath::Abs(1.0 / tpcseed->get_qOverR()));
0196 float beta = std::atan(tpcseed->get_slope());
0197 float alphacorr = std::cos(alpha);
0198 if (alphacorr < 0 || alphacorr > 4)
0199 {
0200 alphacorr = 4;
0201 }
0202 float betacorr = std::cos(beta);
0203 if (betacorr < 0 || betacorr > 4)
0204 {
0205 betacorr = 4;
0206 }
0207 adc /= thick;
0208 adc *= alphacorr;
0209 adc *= betacorr;
0210 dedxlist.push_back(adc);
0211 sort(dedxlist.begin(), dedxlist.end());
0212 }
0213 int trunc_min = 0;
0214 if (dedxlist.empty())
0215 {
0216 return std::numeric_limits<float>::quiet_NaN();
0217 }
0218 int trunc_max = (int) dedxlist.size() * 0.7;
0219 float sumdedx = 0;
0220 int ndedx = 0;
0221 for (int j = trunc_min; j <= trunc_max; j++)
0222 {
0223 sumdedx += dedxlist.at(j);
0224 ndedx++;
0225 }
0226 sumdedx /= ndedx;
0227 return sumdedx;
0228 }
0229
0230 float TpcSeedsQA::cal_track_length(SvtxTrack *track)
0231 {
0232 float minR = std::numeric_limits<float>::max();
0233 float maxR = 0;
0234 for (const auto &ckey : get_cluster_keys(track))
0235 {
0236 auto *cluster = clustermap->findCluster(ckey);
0237
0238
0239 Acts::Vector3 global = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, track->get_crossing());
0240
0241
0242 float R = std::sqrt(pow(global.x(), 2) + pow(global.y(), 2));
0243 minR = std::min(R, minR);
0244 maxR = std::max(R, maxR);
0245 }
0246 float tracklength = maxR - minR;
0247 return tracklength;
0248 }
0249
0250 float *TpcSeedsQA::cal_dedx_cluster(SvtxTrack *track)
0251 {
0252
0253 std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> global_raw;
0254 std::vector<std::pair<TrkrDefs::cluskey, Acts::Vector3>> global_moved;
0255 float minR = std::numeric_limits<float>::max();
0256 float maxR = 0;
0257 for (const auto &ckey : get_cluster_keys(track))
0258 {
0259 auto *cluster = clustermap->findCluster(ckey);
0260
0261
0262 Acts::Vector3 global = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, track->get_crossing());
0263
0264
0265 global_raw.emplace_back(ckey, global);
0266 float R = std::sqrt(pow(global.x(), 2) + pow(global.y(), 2));
0267 minR = std::min(R, minR);
0268 maxR = std::max(R, maxR);
0269 }
0270 float tracklength = maxR - minR;
0271 if (collision_or_cosmics == true && tracklength < 25)
0272 {
0273 float *dedxarray = new float[10];
0274 for (int i = 0; i < 10; ++i)
0275 {
0276 dedxarray[i] = -1;
0277 }
0278 return dedxarray;
0279 }
0280
0281
0282 global_moved = m_clusterMover.processTrack(global_raw);
0283
0284 float fcorr = std::fabs(std::sin(eta_to_theta(track->get_eta())));
0285 Acts::Vector3 clusglob_moved(0, 0, 0);
0286 float adc_z0 = 0;
0287 int nclus_z0 = 0;
0288 float adc_z1 = 0;
0289 int nclus_z1 = 0;
0290 float adc_z2 = 0;
0291 int nclus_z2 = 0;
0292 float adc_z3 = 0;
0293 int nclus_z3 = 0;
0294 float adc_z4 = 0;
0295 int nclus_z4 = 0;
0296 float adc_z5 = 0;
0297 int nclus_z5 = 0;
0298 float adc_z6 = 0;
0299 int nclus_z6 = 0;
0300 float adc_z7 = 0;
0301 int nclus_z7 = 0;
0302 float adc_z8 = 0;
0303 int nclus_z8 = 0;
0304 float adc_z9 = 0;
0305 int nclus_z9 = 0;
0306 for (const auto &pair : global_moved)
0307 {
0308 auto ckey = pair.first;
0309 auto *cluster = clustermap->findCluster(ckey);
0310 clusglob_moved = pair.second;
0311
0312 auto detid = TrkrDefs::getTrkrId(ckey);
0313 if (detid != TrkrDefs::TrkrId::tpcId)
0314 {
0315 continue;
0316 }
0317
0318
0319 auto layer = TrkrDefs::getLayer(ckey);
0320 if (layer < 23)
0321 {
0322 continue;
0323 }
0324
0325 float clusgz = clusglob_moved.z();
0326 if (clusgz > -100 && clusgz <= -80)
0327 {
0328 adc_z0 += cluster->getAdc() * fcorr;
0329 nclus_z0++;
0330 }
0331 else if (clusgz > -80 && clusgz <= -60)
0332 {
0333 adc_z1 += cluster->getAdc() * fcorr;
0334 nclus_z1++;
0335 }
0336 else if (clusgz > -60 && clusgz <= -40)
0337 {
0338 adc_z2 += cluster->getAdc() * fcorr;
0339 nclus_z2++;
0340 }
0341 else if (clusgz > -40 && clusgz <= -20)
0342 {
0343 adc_z3 += cluster->getAdc() * fcorr;
0344 nclus_z3++;
0345 }
0346 else if (clusgz > -20 && clusgz <= 0)
0347 {
0348 adc_z4 += cluster->getAdc() * fcorr;
0349 nclus_z4++;
0350 }
0351 else if (clusgz > 0 && clusgz <= 20)
0352 {
0353 adc_z5 += cluster->getAdc() * fcorr;
0354 nclus_z5++;
0355 }
0356 else if (clusgz > 20 && clusgz <= 40)
0357 {
0358 adc_z6 += cluster->getAdc() * fcorr;
0359 nclus_z6++;
0360 }
0361 else if (clusgz > 40 && clusgz <= 60)
0362 {
0363 adc_z7 += cluster->getAdc() * fcorr;
0364 nclus_z7++;
0365 }
0366 else if (clusgz > 60 && clusgz <= 80)
0367 {
0368 adc_z8 += cluster->getAdc() * fcorr;
0369 nclus_z8++;
0370 }
0371 else if (clusgz > 80 && clusgz <= 100)
0372 {
0373 adc_z9 += cluster->getAdc() * fcorr;
0374 nclus_z9++;
0375 }
0376 }
0377
0378 adc_z0 /= nclus_z0;
0379 adc_z1 /= nclus_z1;
0380 adc_z2 /= nclus_z2;
0381 adc_z3 /= nclus_z3;
0382 adc_z4 /= nclus_z4;
0383 adc_z5 /= nclus_z5;
0384 adc_z6 /= nclus_z6;
0385 adc_z7 /= nclus_z7;
0386 adc_z8 /= nclus_z8;
0387 adc_z9 /= nclus_z9;
0388
0389 float *dedxarray = new float[10];
0390 dedxarray[0] = adc_z0;
0391 dedxarray[1] = adc_z1;
0392 dedxarray[2] = adc_z2;
0393 dedxarray[3] = adc_z3;
0394 dedxarray[4] = adc_z4;
0395 dedxarray[5] = adc_z5;
0396 dedxarray[6] = adc_z6;
0397 dedxarray[7] = adc_z7;
0398 dedxarray[8] = adc_z8;
0399 dedxarray[9] = adc_z9;
0400
0401 return dedxarray;
0402 }
0403
0404
0405 int TpcSeedsQA::process_event(PHCompositeNode *topNode)
0406 {
0407 auto *gl1 = findNode::getClass<Gl1RawHit>(topNode, "GL1RAWHIT");
0408 if (gl1)
0409 {
0410 m_bco = gl1->get_bco();
0411 }
0412 else
0413 {
0414 Gl1Packet *gl1PacketInfo = findNode::getClass<Gl1Packet>(topNode, "GL1RAWHIT");
0415 if (!gl1PacketInfo)
0416 {
0417 m_bco = std::numeric_limits<uint64_t>::quiet_NaN();
0418 }
0419 if (gl1PacketInfo)
0420 {
0421 m_bco = gl1PacketInfo->getBCO();
0422 }
0423 }
0424
0425 h_ntrack1d->Fill(trackmap->size());
0426
0427 std::pair<int, int> ntrack_isfromvtx_pos;
0428 std::pair<int, int> ntrack_isfromvtx_neg;
0429
0430 int ntrack1d_pos = 0;
0431 int ntrack1d_neg = 0;
0432 int ntrack1d_ptg1_pos = 0;
0433 int ntrack1d_ptg1_neg = 0;
0434
0435 int nclus[2][3][12] = {{{0}}};
0436 int madc[2][3][12] = {{{0}}};
0437 for (const auto &[key, track] : *trackmap)
0438 {
0439 if (!track)
0440 {
0441 continue;
0442 }
0443
0444 int charge = track->get_charge();
0445 float quality = track->get_quality();
0446 float pt = track->get_pt();
0447
0448 h_pt->Fill(pt);
0449 if (charge == 1)
0450 {
0451 ntrack1d_pos++;
0452 if (pt > 1)
0453 {
0454 ntrack1d_ptg1_pos++;
0455 }
0456 h_pt_pos->Fill(pt);
0457 }
0458 else if (charge == -1)
0459 {
0460 ntrack1d_neg++;
0461 if (pt > 1)
0462 {
0463 ntrack1d_ptg1_neg++;
0464 }
0465 h_pt_neg->Fill(pt);
0466 }
0467
0468 auto ckeys = get_cluster_keys(track);
0469 std::vector<Acts::Vector3> cluspos;
0470 TrackFitUtils::getTrackletClusters(actsgeom, clustermap, cluspos, ckeys);
0471 float eta = track->get_eta();
0472 float phi = track->get_phi();
0473
0474
0475
0476
0477
0478 int ntpc = 0;
0479 int ntpc_phisize1 = 0;
0480 int nmms = 0;
0481
0482 for (auto &ckey : ckeys)
0483 {
0484 if (TrkrDefs::getTrkrId(ckey) == TrkrDefs::tpcId)
0485 {
0486 TrkrCluster *cluster = clustermap->findCluster(ckey);
0487 if (cluster->getPhiSize() == 1)
0488 {
0489 ntpc_phisize1++;
0490 }
0491 }
0492 switch (TrkrDefs::getTrkrId(ckey))
0493 {
0494
0495
0496
0497
0498
0499
0500 case TrkrDefs::tpcId:
0501 ntpc++;
0502 break;
0503 case TrkrDefs::micromegasId:
0504 nmms++;
0505 break;
0506 default:
0507 break;
0508 }
0509 }
0510
0511 Acts::Vector3 zero = Acts::Vector3::Zero();
0512 auto dcapair_origin = TrackAnalysisUtils::get_dca(track, zero);
0513
0514 auto *trackvtx = vertexmap->get(track->get_vertex_id());
0515 if (!trackvtx)
0516 {
0517
0518 if (charge == 1)
0519 {
0520 ntrack_isfromvtx_pos.first++;
0521 }
0522 else if (charge == -1)
0523 {
0524 ntrack_isfromvtx_neg.first++;
0525 }
0526 }
0527 else
0528 {
0529 Acts::Vector3 track_vtx(trackvtx->get_x(), trackvtx->get_y(), trackvtx->get_z());
0530 auto dcapair_vtx = TrackAnalysisUtils::get_dca(track, track_vtx);
0531 if (charge == 1)
0532 {
0533 ntrack_isfromvtx_pos.second++;
0534 h_dcaxyvtx_phi_pos->Fill(phi, dcapair_vtx.first.first);
0535 h_dcazvtx_phi_pos->Fill(phi, dcapair_vtx.second.first);
0536 }
0537 else if (charge == -1)
0538 {
0539 ntrack_isfromvtx_neg.second++;
0540 h_dcaxyvtx_phi_neg->Fill(phi, dcapair_vtx.first.first);
0541 h_dcazvtx_phi_neg->Fill(phi, dcapair_vtx.second.first);
0542 }
0543 }
0544
0545 if (charge == 1)
0546 {
0547 h_ntpc_fullpt_pos->Fill(ntpc);
0548 if (dcapair_origin.second.first > 0)
0549 {
0550 h_dcaxyorigin_phi_north_pos->Fill(phi, dcapair_origin.first.first);
0551 }
0552 else if (dcapair_origin.second.first <= 0)
0553 {
0554 h_dcaxyorigin_phi_south_pos->Fill(phi, dcapair_origin.first.first);
0555 }
0556 h_dcazorigin_phi_pos->Fill(phi, dcapair_origin.second.first);
0557 if (pt > 1)
0558 {
0559 h_ntrack_pos->Fill(eta, phi);
0560 if (trackvtx)
0561 {
0562 float vz = trackvtx->get_z();
0563 float eta_min = cal_tpc_eta_min_max(vz).first;
0564 float eta_max = cal_tpc_eta_min_max(vz).second;
0565 if (eta > eta_min && eta < eta_max)
0566 {
0567 h_ntpc_pos->Fill(ntpc);
0568 }
0569 }
0570 h_ntpot_pos->Fill(nmms);
0571 h_ntpc_quality_pos->Fill(ntpc, quality);
0572 h_avgnclus_eta_phi_pos->Fill(eta, phi, ntpc);
0573
0574 h_cluster_phisize1_fraction_pos->Fill((double) ntpc_phisize1 / (double) ntpc);
0575 }
0576 }
0577 else if (charge == -1)
0578 {
0579 h_ntpc_fullpt_neg->Fill(ntpc);
0580 if (dcapair_origin.second.first > 0)
0581 {
0582 h_dcaxyorigin_phi_north_neg->Fill(phi, dcapair_origin.first.first);
0583 }
0584 else if (dcapair_origin.second.first <= 0)
0585 {
0586 h_dcaxyorigin_phi_south_neg->Fill(phi, dcapair_origin.first.first);
0587 }
0588 h_dcazorigin_phi_neg->Fill(phi, dcapair_origin.second.first);
0589 if (pt > 1)
0590 {
0591 h_ntrack_neg->Fill(eta, phi);
0592 if (trackvtx)
0593 {
0594 float vz = trackvtx->get_z();
0595 float eta_min = cal_tpc_eta_min_max(vz).first;
0596 float eta_max = cal_tpc_eta_min_max(vz).second;
0597 if (eta > eta_min && eta < eta_max)
0598 {
0599 h_ntpc_neg->Fill(ntpc);
0600 }
0601 }
0602 h_ntpot_neg->Fill(nmms);
0603 h_ntpc_quality_neg->Fill(ntpc, quality);
0604 h_avgnclus_eta_phi_neg->Fill(eta, phi, ntpc);
0605
0606 h_cluster_phisize1_fraction_neg->Fill((double) ntpc_phisize1 / (double) ntpc);
0607 }
0608 }
0609 }
0610 h_ntrack1d_pos->Fill(ntrack1d_pos);
0611 h_ntrack1d_neg->Fill(ntrack1d_neg);
0612 h_ntrack1d_ptg1_pos->Fill(ntrack1d_ptg1_pos);
0613 h_ntrack1d_ptg1_neg->Fill(ntrack1d_ptg1_neg);
0614 h_ntrack1d_ptg1->Fill(ntrack1d_ptg1_pos + ntrack1d_ptg1_neg);
0615
0616 h_ntrack_isfromvtx_pos->SetBinContent(1, h_ntrack_isfromvtx_pos->GetBinContent(1) + ntrack_isfromvtx_pos.first);
0617 h_ntrack_isfromvtx_pos->SetBinContent(2, h_ntrack_isfromvtx_pos->GetBinContent(2) + ntrack_isfromvtx_pos.second);
0618 h_ntrack_isfromvtx_neg->SetBinContent(1, h_ntrack_isfromvtx_neg->GetBinContent(1) + ntrack_isfromvtx_neg.first);
0619 h_ntrack_isfromvtx_neg->SetBinContent(2, h_ntrack_isfromvtx_neg->GetBinContent(2) + ntrack_isfromvtx_neg.second);
0620
0621
0622 h_nvertex->Fill(vertexmap->size());
0623 for (const auto &[key, vertex] : *vertexmap)
0624 {
0625 if (!vertex)
0626 {
0627 continue;
0628 }
0629
0630 float vx = vertex->get_x();
0631 float vy = vertex->get_y();
0632 float vz = vertex->get_z();
0633 float vt = vertex->get_t0();
0634 float vchi2 = vertex->get_chisq();
0635 int vndof = vertex->get_ndof();
0636
0637
0638
0639
0640 h_vx->Fill(vx);
0641 h_vy->Fill(vy);
0642 h_vx_vy->Fill(vx, vy);
0643 h_vz->Fill(vz);
0644 h_vt->Fill(vt);
0645 h_vchi2dof->Fill(float(vchi2 / vndof));
0646
0647
0648 h_ntrackpervertex->Fill(vertex->size_tracks());
0649 }
0650
0651 auto fill = [](TH1 *h, float val)
0652 { if (h) { h->Fill(val); } };
0653
0654 std::set<unsigned int> tpc_seed_ids;
0655 for (const auto &[key, track] : *trackmap)
0656 {
0657 if (!track)
0658 {
0659 continue;
0660 }
0661 m_px = track->get_px();
0662 m_py = track->get_py();
0663 m_pz = track->get_pz();
0664 m_pt = std::sqrt((m_px * m_px) + (m_py * m_py));
0665 m_ptot = std::sqrt((m_px * m_px) + (m_py * m_py) + (m_pz * m_pz));
0666 TrackSeed *tpcseed = track->get_tpc_seed();
0667 m_charge = track->get_charge();
0668 m_dedx = calc_dedx(tpcseed);
0669
0670 m_ntpc = 0;
0671 m_region.clear();
0672 m_clusgz.clear();
0673 m_cluslayer.clear();
0674 m_clusphisize.clear();
0675 m_cluszsize.clear();
0676
0677 for (const auto &ckey : get_cluster_keys(track))
0678 {
0679 if (TrkrDefs::getTrkrId(ckey) == TrkrDefs::tpcId)
0680 {
0681 m_ntpc++;
0682 }
0683 }
0684
0685 for (const auto &ckey : get_cluster_keys(track))
0686 {
0687 TrkrCluster *cluster = clustermap->findCluster(ckey);
0688 const Acts::Vector3 clusglob = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, track->get_crossing());
0689
0690 switch (TrkrDefs::getTrkrId(ckey))
0691 {
0692 case TrkrDefs::tpcId:
0693 const auto it = m_layerRegionMap.find(TrkrDefs::getLayer(ckey));
0694 int region = it->second;
0695 m_region.push_back(region);
0696 m_clusgz.push_back(clusglob.z());
0697 m_cluslayer.push_back(TrkrDefs::getLayer(ckey));
0698 m_clusphisize.push_back(cluster->getPhiSize());
0699 m_cluszsize.push_back(cluster->getZSize());
0700 int this_sector = (int) TpcDefs::getSectorId(ckey);
0701 int this_side = (int) TpcDefs::getSide(ckey);
0702 int is_onepad = 0;
0703 if (cluster->getPhiSize() <= 1)
0704 {
0705 is_onepad = 1;
0706 }
0707 if (m_ntpc > 30 && cluster->getPhiSize() > 1 && m_dedx < 1500 && m_pt > 1.0 && m_pt < 50)
0708 {
0709 h_adc_sector[region]->Fill((this_sector + 1) * (2 * (this_side - 0.5)), cluster->getAdc());
0710 }
0711 if (m_ntpc > 30 && m_dedx < 1000 && m_pt > 1.0 && m_pt < 50 && m_charge < 0)
0712 {
0713 h_onepad_frac[region]->Fill((this_sector + 1) * (2 * (this_side - 0.5)), is_onepad);
0714 }
0715 nclus[this_side][region][this_sector] += 1;
0716 madc[this_side][region][this_sector] += cluster->getAdc();
0717 break;
0718 }
0719 }
0720
0721 if (m_ptot > 0.2 && m_ptot < 4 && m_ntpc > 30)
0722 {
0723 h_dedx->Fill(m_charge * m_ptot, m_dedx);
0724 if (collision_or_cosmics == false || (collision_or_cosmics == true && cal_track_length(track) > 25))
0725 {
0726 h_dedx_pcaz->Fill(track->get_z(), m_dedx);
0727 }
0728 }
0729
0730 if (m_ptot > 1.0 && m_pt < 4 && m_ntpc > 30 && m_charge < 0 && m_dedx < 1000 && m_dedx > 50)
0731 {
0732 h_mip_dedx->Fill(m_dedx);
0733 }
0734
0735 if (m_ntpc > 30)
0736 {
0737 float *cluster_dedx = cal_dedx_cluster(track);
0738 for (int iz = 0; iz < 10; iz++)
0739 {
0740 h_dedx_pq_z[iz]->Fill(m_charge * m_ptot, cluster_dedx[iz]);
0741 }
0742 }
0743
0744
0745
0746
0747
0748
0749 int nClus = m_cluslayer.size();
0750 for (int cl = 0; cl < nClus; cl++)
0751 {
0752 if (m_pt > 1 && m_ntpc > 25)
0753 {
0754 if (m_clusphisize[cl] == 1 && m_cluszsize[cl] > 1)
0755 {
0756 if (m_clusgz[cl] < 0.)
0757 {
0758 const auto hiter = phihistos.find(m_region[cl]);
0759 if (hiter == phihistos.end())
0760 {
0761 continue;
0762 }
0763 fill(hiter->second.cphisize1pT_side0, m_pt);
0764 }
0765 else if (m_clusgz[cl] > 0.)
0766 {
0767 const auto hiter = phihistos.find(m_region[cl]);
0768 if (hiter == phihistos.end())
0769 {
0770 continue;
0771 }
0772 fill(hiter->second.cphisize1pT_side1, m_pt);
0773 }
0774 }
0775 if (m_clusphisize[cl] >= 1 && m_cluszsize[cl] > 1)
0776 {
0777 if (m_clusgz[cl] < 0.)
0778 {
0779 const auto hiter = phihistos.find(m_region[cl]);
0780 if (hiter == phihistos.end())
0781 {
0782 continue;
0783 }
0784 fill(hiter->second.cphisizegeq1pT_side0, m_pt);
0785 }
0786 else if (m_clusgz[cl] > 0.)
0787 {
0788 const auto hiter = phihistos.find(m_region[cl]);
0789 if (hiter == phihistos.end())
0790 {
0791 continue;
0792 }
0793 fill(hiter->second.cphisizegeq1pT_side1, m_pt);
0794 }
0795 }
0796 }
0797 }
0798
0799 for (auto &pair : phihistos)
0800 {
0801 pair.second.Clear();
0802 }
0803
0804 for (int cl = 0; cl < nClus; cl++)
0805 {
0806 if (m_pt > 1 && m_ntpc > 25)
0807 {
0808 if (m_clusgz[cl] < 0.)
0809 {
0810 const auto hiter = phihistos.find(m_region[cl]);
0811 if (hiter == phihistos.end())
0812 {
0813 continue;
0814 }
0815 if (m_clusphisize[cl] == 1 && m_cluszsize[cl] > 1)
0816 {
0817 hiter->second.ntpc_side0_phisize1++;
0818 }
0819 if (m_clusphisize[cl] >= 1 && m_cluszsize[cl] > 1)
0820 {
0821 hiter->second.ntpc_side0++;
0822 }
0823 }
0824 else if (m_clusgz[cl] > 0.)
0825 {
0826 const auto hiter = phihistos.find(m_region[cl]);
0827 if (hiter == phihistos.end())
0828 {
0829 continue;
0830 }
0831 if (m_clusphisize[cl] == 1 && m_cluszsize[cl] > 1)
0832 {
0833 hiter->second.ntpc_side1_phisize1++;
0834 }
0835 if (m_clusphisize[cl] >= 1 && m_cluszsize[cl] > 1)
0836 {
0837 hiter->second.ntpc_side1++;
0838 }
0839 }
0840 }
0841 }
0842
0843 for (const auto ®ion : {0, 1, 2})
0844 {
0845 if (phihistos[region].ntpc_side0 > 0)
0846 {
0847 double frac_side0 = (double) phihistos[region].ntpc_side0_phisize1 / (double) phihistos[region].ntpc_side0;
0848 h_cluster_phisize1_fraction_side0[region]->Fill(frac_side0);
0849 h_cluster_phisize1_fraction_pt_side0[region]->Fill(m_pt, frac_side0);
0850
0851 int index_pt_side0 = h_cluster_phisize1_fraction_mean_side0[region]->FindBin(m_pt) - 1;
0852 if (index_pt_side0 < h_cluster_phisize1_fraction_mean_side0[region]->GetNbinsX())
0853 {
0854 frac_side0_pt[region][index_pt_side0] += frac_side0;
0855 num_track_side0_pt[region][index_pt_side0]++;
0856 }
0857 }
0858
0859 if (phihistos[region].ntpc_side1 > 0)
0860 {
0861 double frac_side1 = (double) phihistos[region].ntpc_side1_phisize1 / (double) phihistos[region].ntpc_side1;
0862 h_cluster_phisize1_fraction_side1[region]->Fill(frac_side1);
0863 h_cluster_phisize1_fraction_pt_side1[region]->Fill(m_pt, frac_side1);
0864
0865 int index_pt_side1 = h_cluster_phisize1_fraction_mean_side1[region]->FindBin(m_pt) - 1;
0866 if (index_pt_side1 < h_cluster_phisize1_fraction_mean_side1[region]->GetNbinsX())
0867 {
0868 frac_side1_pt[region][index_pt_side1] += frac_side1;
0869 num_track_side1_pt[region][index_pt_side1]++;
0870 }
0871 }
0872 }
0873 }
0874
0875 for (int iside = 0; iside < 2; iside++)
0876 {
0877 for (int iregion = 0; iregion < 3; iregion++)
0878 {
0879 for (int isector = 0; isector < 12; isector++)
0880 {
0881
0882 int nCluster = nclus[iside][iregion][isector];
0883 float meanAdc = 0;
0884 if (nCluster > 0)
0885 {
0886 meanAdc = madc[iside][iregion][isector] / (nCluster * 1.);
0887 }
0888 nt_sector_event_summary->Fill(m_event, m_segment, m_bco, iside, iregion, isector, nCluster, meanAdc);
0889 }
0890 }
0891 }
0892 m_event++;
0893 return Fun4AllReturnCodes::EVENT_OK;
0894 }
0895
0896 std::vector<TrkrDefs::cluskey> TpcSeedsQA::get_cluster_keys(SvtxTrack *track)
0897 {
0898 std::vector<TrkrDefs::cluskey> out;
0899 for (const auto &seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0900 {
0901 if (seed)
0902 {
0903 std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0904 }
0905 }
0906 return out;
0907 }
0908
0909
0910 int TpcSeedsQA::EndRun(const int )
0911 {
0912 for (const auto ®ion : {0, 1, 2})
0913 {
0914 for (const auto &index_pt : {0, 1, 2, 3})
0915 {
0916 if (num_track_side0_pt[region][index_pt] > 0)
0917 {
0918 h_cluster_phisize1_fraction_mean_side0[region]->SetBinContent(index_pt + 1, frac_side0_pt[region][index_pt] / num_track_side0_pt[region][index_pt]);
0919 h_cluster_phisize1_fraction_mean_numerator_side0[region]->SetBinContent(index_pt + 1, frac_side0_pt[region][index_pt]);
0920 h_cluster_phisize1_fraction_mean_denominator_side0[region]->SetBinContent(index_pt + 1, num_track_side0_pt[region][index_pt]);
0921 }
0922 if (num_track_side1_pt[region][index_pt] > 0)
0923 {
0924 h_cluster_phisize1_fraction_mean_side1[region]->SetBinContent(index_pt + 1, frac_side1_pt[region][index_pt] / num_track_side1_pt[region][index_pt]);
0925 h_cluster_phisize1_fraction_mean_numerator_side1[region]->SetBinContent(index_pt + 1, frac_side1_pt[region][index_pt]);
0926 h_cluster_phisize1_fraction_mean_denominator_side1[region]->SetBinContent(index_pt + 1, num_track_side1_pt[region][index_pt]);
0927 }
0928 }
0929 }
0930
0931 auto *hm = QAHistManagerDef::getHistoManager();
0932 assert(hm);
0933
0934 return Fun4AllReturnCodes::EVENT_OK;
0935 }
0936
0937
0938 int TpcSeedsQA::End(PHCompositeNode * )
0939 {
0940 return Fun4AllReturnCodes::EVENT_OK;
0941 }
0942
0943 std::string TpcSeedsQA::getHistoPrefix() const
0944 {
0945 return std::string("h_") + Name() + std::string("_");
0946 }
0947
0948 void TpcSeedsQA::createHistos()
0949 {
0950 auto *hm = QAHistManagerDef::getHistoManager();
0951 assert(hm);
0952
0953 {
0954 auto *h = new TH1F(std::string(getHistoPrefix() + "ntpc_fullpt_pos").c_str(), "TPC clusters per positive track;Number of TPC clusters per positive track;Entries", 55, -0.5, 54.5);
0955 hm->registerHisto(h);
0956 }
0957
0958 {
0959 auto *h = new TH1F(std::string(getHistoPrefix() + "ntpc_fullpt_neg").c_str(), "TPC clusters per negative track;Number of TPC clusters per negative track;Entries", 55, -0.5, 54.5);
0960 hm->registerHisto(h);
0961 }
0962
0963 {
0964 auto *h = new TH1F(std::string(getHistoPrefix() + "ntpc_pos").c_str(), "TPC clusters per positive track (pT>1GeV,eta cut);Number of TPC clusters per positive track;Entries", 55, -0.5, 54.5);
0965 hm->registerHisto(h);
0966 }
0967
0968 {
0969 auto *h = new TH1F(std::string(getHistoPrefix() + "ntpc_neg").c_str(), "TPC clusters per negative track (pT>1GeV,eta cut);Number of TPC clusters per negative track;Entries", 55, -0.5, 54.5);
0970 hm->registerHisto(h);
0971 }
0972
0973 {
0974 auto *h = new TH1F(std::string(getHistoPrefix() + "ntpot_pos").c_str(), "TPOT clusters per positive track (pT>1GeV);Number of TPOT clusters per positive track;Entries", 2, -0.5, 1.5);
0975 hm->registerHisto(h);
0976 }
0977
0978 {
0979 auto *h = new TH1F(std::string(getHistoPrefix() + "ntpot_neg").c_str(), "TPOT clusters per negative track (pT>1GeV);Number of TPOT clusters per negative track;Entries", 2, -0.5, 1.5);
0980 hm->registerHisto(h);
0981 }
0982
0983 {
0984 auto *h = new TH2F(std::string(getHistoPrefix() + "ntpc_quality_pos").c_str(), "Number of TPC clusters per positive track (pT>1GeV);Number of TPC clusters per positive track;Quality", 55, -0.5, 54.5, 100, 0, 10);
0985 hm->registerHisto(h);
0986 }
0987
0988 {
0989 auto *h = new TH2F(std::string(getHistoPrefix() + "ntpc_quality_neg").c_str(), "Number of TPC clusters per negative track (pT>1GeV);Number of TPC clusters per negative track;Quality", 55, -0.5, 54.5, 100, 0, 10);
0990 hm->registerHisto(h);
0991 }
0992
0993 {
0994 auto *h = new TH1F(std::string(getHistoPrefix() + "nrecotracks1d").c_str(), "Number of reconstructed tracks;Number of TPC tracklets;Entries", 50, 0, 200);
0995 hm->registerHisto(h);
0996 }
0997
0998 {
0999 auto *h = new TH1F(std::string(getHistoPrefix() + "nrecotracks1d_pos").c_str(), "Number of reconstructed positive tracks;Number of positive TPC tracklets;Entries", 50, 0, 200);
1000 hm->registerHisto(h);
1001 }
1002
1003 {
1004 auto *h = new TH1F(std::string(getHistoPrefix() + "nrecotracks1d_neg").c_str(), "Number of reconstructed negative tracks;Number of negative TPC tracklets;Entries", 50, 0, 200);
1005 hm->registerHisto(h);
1006 }
1007
1008 {
1009 auto *h = new TH1F(std::string(getHistoPrefix() + "nrecotracks1d_ptg1").c_str(), "Number of reconstructed tracks (pT>1GeV);Number of TPC tracklets;Entries", 50, 0, 200);
1010 hm->registerHisto(h);
1011 }
1012
1013 {
1014 auto *h = new TH1F(std::string(getHistoPrefix() + "nrecotracks1d_ptg1_pos").c_str(), "Number of reconstructed positive tracks (pT>1GeV);Number of positive TPC tracklets;Entries", 50, 0, 200);
1015 hm->registerHisto(h);
1016 }
1017
1018 {
1019 auto *h = new TH1F(std::string(getHistoPrefix() + "nrecotracks1d_ptg1_neg").c_str(), "Number of reconstructed negative tracks (pT>1GeV);Number of negative TPC tracklets;Entries", 50, 0, 200);
1020 hm->registerHisto(h);
1021 }
1022
1023 {
1024 auto *h = new TH1F(std::string(getHistoPrefix() + "pt").c_str(), "p_{T} distribution of reconstructed tracks;Track p_{T};Entries", 100, 0, 10);
1025 hm->registerHisto(h);
1026 }
1027
1028 {
1029 auto *h = new TH1F(std::string(getHistoPrefix() + "pt_pos").c_str(), "p_{T} distribution of reconstructed positive tracks;Track p_{T};Entries", 100, 0, 10);
1030 hm->registerHisto(h);
1031 }
1032
1033 {
1034 auto *h = new TH1F(std::string(getHistoPrefix() + "pt_neg").c_str(), "p_{T} distribution of reconstructed negative tracks;Track p_{T};Entries", 100, 0, 10);
1035 hm->registerHisto(h);
1036 }
1037
1038 {
1039 auto *h = new TH2F(std::string(getHistoPrefix() + "nrecotracks_pos").c_str(), "Number of reconstructed positive tracks (pT>1GeV);#eta;#phi [rad];Entries", 100, -1.1, 1.1, 300, -3.14159, 3.1459);
1040 hm->registerHisto(h);
1041 }
1042
1043 {
1044 auto *h = new TH2F(std::string(getHistoPrefix() + "nrecotracks_neg").c_str(), "Number of reconstructed negative tracks (pT>1GeV);#eta;#phi [rad];Entries", 100, -1.1, 1.1, 300, -3.14159, 3.1459);
1045 hm->registerHisto(h);
1046 }
1047
1048 {
1049 auto *h = new TProfile2D(std::string(getHistoPrefix() + "avgnclus_eta_phi_pos").c_str(), "Average number of clusters per positive track (pT>1GeV);#eta;#phi [rad];Average number of clusters per positive track", 100, -1.1, 1.1, 300, -3.14159, 3.1459, 0, 55);
1050 hm->registerHisto(h);
1051 }
1052
1053 {
1054 auto *h = new TProfile2D(std::string(getHistoPrefix() + "avgnclus_eta_phi_neg").c_str(), "Average number of clusters per negative track (pT>1GeV);#eta;#phi [rad];Average number of clusters per negative track", 100, -1.1, 1.1, 300, -3.14159, 3.1459, 0, 55);
1055 hm->registerHisto(h);
1056 }
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068 {
1069 auto *h = new TH2F(std::string(getHistoPrefix() + "dcaxyorigin_phi_north_pos").c_str(), "DCA xy origin vs phi for positive track (dcaz>0);#phi [rad];DCA_{xy} wrt origin [cm];Entries", 300, -3.14159, 3.1459, 100, -10, 10);
1070 hm->registerHisto(h);
1071 }
1072
1073 {
1074 auto *h = new TH2F(std::string(getHistoPrefix() + "dcaxyorigin_phi_south_pos").c_str(), "DCA xy origin vs phi for positive track (dcaz<0);#phi [rad];DCA_{xy} wrt origin [cm];Entries", 300, -3.14159, 3.1459, 100, -10, 10);
1075 hm->registerHisto(h);
1076 }
1077
1078 {
1079 auto *h = new TH2F(std::string(getHistoPrefix() + "dcaxyorigin_phi_north_neg").c_str(), "DCA xy origin vs phi for negative track (dcaz>0);#phi [rad];DCA_{xy} wrt origin [cm];Entries", 300, -3.14159, 3.1459, 100, -10, 10);
1080 hm->registerHisto(h);
1081 }
1082
1083 {
1084 auto *h = new TH2F(std::string(getHistoPrefix() + "dcaxyorigin_phi_south_neg").c_str(), "DCA xy origin vs phi for negative track (dcaz<0);#phi [rad];DCA_{xy} wrt origin [cm];Entries", 300, -3.14159, 3.1459, 100, -10, 10);
1085 hm->registerHisto(h);
1086 }
1087
1088 {
1089 auto *h = new TH2F(std::string(getHistoPrefix() + "dcaxyvtx_phi_pos").c_str(), "DCA xy vertex vs phi for positive track;#phi [rad];DCA_{xy} wrt vertex [cm];Entries", 300, -3.14159, 3.1459, 90, -3, 3);
1090 hm->registerHisto(h);
1091 }
1092
1093 {
1094 auto *h = new TH2F(std::string(getHistoPrefix() + "dcaxyvtx_phi_neg").c_str(), "DCA xy vertex vs phi for negative track;#phi [rad];DCA_{xy} wrt vertex [cm];Entries", 300, -3.14159, 3.1459, 90, -3, 3);
1095 hm->registerHisto(h);
1096 }
1097
1098 {
1099 auto *h = new TH2F(std::string(getHistoPrefix() + "dcazorigin_phi_pos").c_str(), "DCA z origin vs phi for positive track;#phi [rad];DCA_{z} wrt origin [cm];Entries", 300, -3.14159, 3.1459, 100, -10, 10);
1100 hm->registerHisto(h);
1101 }
1102
1103 {
1104 auto *h = new TH2F(std::string(getHistoPrefix() + "dcazorigin_phi_neg").c_str(), "DCA z origin vs phi for negative track;#phi [rad];DCA_{z} wrt origin [cm];Entries", 300, -3.14159, 3.1459, 100, -10, 10);
1105 hm->registerHisto(h);
1106 }
1107
1108 {
1109 auto *h = new TH2F(std::string(getHistoPrefix() + "dcazvtx_phi_pos").c_str(), "DCA z vertex vs phi for positive track;#phi [rad];DCA_{z} wrt vertex [cm];Entries", 300, -3.14159, 3.1459, 100, -10, 10);
1110 hm->registerHisto(h);
1111 }
1112
1113 {
1114 auto *h = new TH2F(std::string(getHistoPrefix() + "dcazvtx_phi_neg").c_str(), "DCA z vertex vs phi for negative track;#phi [rad];DCA_{z} wrt vertex [cm];Entries", 300, -3.14159, 3.1459, 100, -10, 10);
1115 hm->registerHisto(h);
1116 }
1117
1118 {
1119 auto *h = new TH1F(std::string(getHistoPrefix() + "ntrack_isfromvtx_pos").c_str(), "Num of positive tracks associated to a vertex;Is track associated to a vertex;Entries", 2, -0.5, 1.5);
1120 hm->registerHisto(h);
1121 }
1122
1123 {
1124 auto *h = new TH1F(std::string(getHistoPrefix() + "ntrack_isfromvtx_neg").c_str(), "Num of negative tracks associated to a vertex;Is track associated to a vertex;Entries", 2, -0.5, 1.5);
1125 hm->registerHisto(h);
1126 }
1127
1128 {
1129 auto *h = new TH1F(std::string(getHistoPrefix() + "cluster_phisize1_fraction_pos").c_str(), "Fraction of TPC clusters per positive track with phi size of 1 (pT>1GeV);Fraction of TPC clusters phi size of 1;Entries", 100, 0, 1);
1130 hm->registerHisto(h);
1131 }
1132
1133 {
1134 auto *h = new TH1F(std::string(getHistoPrefix() + "cluster_phisize1_fraction_neg").c_str(), "Fraction of TPC clusters per negative track with phi size of 1 (pT>1GeV);Fraction of TPC clusters phi size of 1;Entries", 100, 0, 1);
1135 hm->registerHisto(h);
1136 }
1137
1138
1139 {
1140 auto *h = new TH1F(std::string(getHistoPrefix() + "nrecovertices").c_str(), "Num of reco vertices per event;Number of vertices;Entries", 20, 0, 20);
1141 hm->registerHisto(h);
1142 }
1143
1144 {
1145 auto *h = new TH1F(std::string(getHistoPrefix() + "vx").c_str(), "Vertex x;Vertex x [cm];Entries", 100, -2.5, 2.5);
1146 hm->registerHisto(h);
1147 }
1148
1149 {
1150 auto *h = new TH1F(std::string(getHistoPrefix() + "vy").c_str(), "Vertex y;Vertex y [cm];Entries", 100, -2.5, 2.5);
1151 hm->registerHisto(h);
1152 }
1153
1154 {
1155 auto *h = 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);
1156 hm->registerHisto(h);
1157 }
1158
1159 {
1160 auto *h = new TH1F(std::string(getHistoPrefix() + "vz").c_str(), "Vertex z;Vertex z [cm];Entries", 50, -25, 25);
1161 hm->registerHisto(h);
1162 }
1163
1164 {
1165 auto *h = new TH1F(std::string(getHistoPrefix() + "vt").c_str(), "Vertex t;Vertex t [ns];Entries", 100, -1000, 20000);
1166 hm->registerHisto(h);
1167 }
1168
1169
1170
1171
1172
1173
1174 {
1175 auto *h = new TH1F(std::string(getHistoPrefix() + "vertexchi2dof").c_str(), "Vertex chi2/ndof;Vertex #chi2/ndof;Entries", 100, 0, 20);
1176 hm->registerHisto(h);
1177 }
1178
1179 {
1180 auto *h = new TH1F(std::string(getHistoPrefix() + "ntrackspervertex").c_str(), "Num of tracks per vertex;Number of tracks per vertex;Entries", 50, 0, 50);
1181 hm->registerHisto(h);
1182 }
1183 {
1184 auto *h = new TH2F(std::string(getHistoPrefix() + "dedx").c_str(),
1185 "Num of tracks per vertex;Number of tracks per vertex;Entries",
1186 500, -2, 2, 500, 0, 3000);
1187 hm->registerHisto(h);
1188 }
1189
1190 {
1191 auto *h = new TH1F(std::string(getHistoPrefix() + "mip_dedx").c_str(),
1192 "dEdx of MIPs",
1193 100, 0, 1000);
1194 hm->registerHisto(h);
1195 }
1196
1197 {
1198 auto *h = new TH2F(std::string(getHistoPrefix() + "dedx_pcaz").c_str(),
1199 "track dEdx vs pcaz",
1200 100, -100, 100, 50, 0, 2000);
1201 hm->registerHisto(h);
1202 }
1203
1204 for (int i = 0; i < 10; i++)
1205 {
1206 h_dedx_pq_z[i] = new TH2F(std::format("{}dedx_pq_{}", getHistoPrefix(), i).c_str(),
1207 std::format("mean cluster dEdx at R2 & R3 corrected by path length according to track eta in z-region {}", i).c_str(), 100, -3, 3, 100, 0, 10000);
1208 hm->registerHisto(h_dedx_pq_z[i]);
1209 }
1210
1211 {
1212 auto *nt = new TNtuple(std::string(getHistoPrefix() + "sector_event_summary").c_str(),
1213 "sector_event_summary", "event:segment:bco:side:region:sector:ncluster:meanadc");
1214 hm->registerHisto(nt);
1215 }
1216
1217 for (const auto ®ion : {0, 1, 2})
1218 {
1219 h_adc_sector[region] = new TH2F(std::format("{}adc_sector_{}", getHistoPrefix(), region).c_str(),
1220 std::format("ADC spectrum per, region_{}", region).c_str(), 25, -12.5, 12.5, 50, 0, 1500);
1221 hm->registerHisto(h_adc_sector[region]);
1222
1223 h_onepad_frac[region] = new TProfile(std::format("{}onepad_frac_{}", getHistoPrefix(), region).c_str(),
1224 std::format("TPC Cluster Phi Size == 1 fraction per sector, region_{}", region).c_str(), 25, -12.5, 12.5);
1225 hm->registerHisto(h_onepad_frac[region]);
1226
1227 h_clusphisize1pt_side0[region] = new TH1F(std::format("{}clusphisize1pT_side0_{}", getHistoPrefix(), region).c_str(),
1228 std::format("TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 4, 1, 3.2);
1229 h_clusphisize1pt_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1230 hm->registerHisto(h_clusphisize1pt_side0[region]);
1231
1232 h_clusphisize1pt_side1[region] = new TH1F(std::format("{}clusphisize1pT_side1_{}", getHistoPrefix(), region).c_str(),
1233 std::format("TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 4, 1, 3.2);
1234 h_clusphisize1pt_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1235 hm->registerHisto(h_clusphisize1pt_side1[region]);
1236
1237 h_clusphisizegeq1pt_side0[region] = new TH1F(std::format("{}clusphisizegeq1pT_side0_{}", getHistoPrefix(), region).c_str(),
1238 std::format("TPC Cluster Phi Size >= 1, side 0, region_{}", region).c_str(), 4, 1, 3.2);
1239 h_clusphisizegeq1pt_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1240 hm->registerHisto(h_clusphisizegeq1pt_side0[region]);
1241
1242 h_clusphisizegeq1pt_side1[region] = new TH1F(std::format("{}clusphisizegeq1pT_side1_{}", getHistoPrefix(), region).c_str(),
1243 std::format("TPC Cluster Phi Size >= 1, side 1, region_{}", region).c_str(), 4, 1, 3.2);
1244 h_clusphisizegeq1pt_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1245 hm->registerHisto(h_clusphisizegeq1pt_side1[region]);
1246
1247 h_cluster_phisize1_fraction_side0[region] = new TH1F(std::format("{}sclusphisize1frac_side0_{}", getHistoPrefix(), region).c_str(),
1248 std::format("Fraction of TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 100, 0, 1);
1249 h_cluster_phisize1_fraction_side0[region]->GetXaxis()->SetTitle("Fraction");
1250 hm->registerHisto(h_cluster_phisize1_fraction_side0[region]);
1251
1252 h_cluster_phisize1_fraction_side1[region] = new TH1F(std::format("{}clusphisize1frac_side1_{}", getHistoPrefix(), region).c_str(),
1253 std::format("Fraction of TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 100, 0, 1);
1254 h_cluster_phisize1_fraction_side1[region]->GetXaxis()->SetTitle("Fraction");
1255 hm->registerHisto(h_cluster_phisize1_fraction_side1[region]);
1256
1257 h_cluster_phisize1_fraction_pt_side0[region] = new TH2F(std::format("{}clusphisize1frac_pt_side0_{}", getHistoPrefix(), region).c_str(),
1258 std::format("Pt vs. Fraction of TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 4, 1, 3.2, 100, 0, 1);
1259 h_cluster_phisize1_fraction_pt_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1260 h_cluster_phisize1_fraction_pt_side0[region]->GetYaxis()->SetTitle("Fraction");
1261 hm->registerHisto(h_cluster_phisize1_fraction_pt_side0[region]);
1262
1263 h_cluster_phisize1_fraction_pt_side1[region] = new TH2F(std::format("{}clusphisize1frac_pt_side1_{}", getHistoPrefix(), region).c_str(),
1264 std::format("Pt vs. Fraction of TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 4, 1, 3.2, 100, 0, 1);
1265 h_cluster_phisize1_fraction_pt_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1266 h_cluster_phisize1_fraction_pt_side1[region]->GetYaxis()->SetTitle("Fraction");
1267 hm->registerHisto(h_cluster_phisize1_fraction_pt_side1[region]);
1268
1269 h_cluster_phisize1_fraction_mean_side0[region] = new TH1F(std::format("{}clusphisize1frac_mean_side0_{}", getHistoPrefix(), region).c_str(),
1270 std::format("Pt vs. Average fraction of TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 4, 1, 3.2);
1271 h_cluster_phisize1_fraction_mean_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1272 h_cluster_phisize1_fraction_mean_side0[region]->GetYaxis()->SetTitle("Fraction");
1273 hm->registerHisto(h_cluster_phisize1_fraction_mean_side0[region]);
1274
1275 h_cluster_phisize1_fraction_mean_side1[region] = new TH1F(std::format("{}clusphisize1frac_mean_side1_{}", getHistoPrefix(), region).c_str(),
1276 std::format("Pt vs. Average fraction of TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 4, 1, 3.2);
1277 h_cluster_phisize1_fraction_mean_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1278 h_cluster_phisize1_fraction_mean_side1[region]->GetYaxis()->SetTitle("Fraction");
1279 hm->registerHisto(h_cluster_phisize1_fraction_mean_side1[region]);
1280
1281 h_cluster_phisize1_fraction_mean_numerator_side0[region] = new TH1F(std::format("{}clusphisize1frac_mean_numerator_side0_{}", getHistoPrefix(), region).c_str(),
1282 std::format("Pt vs. Average fraction mean_numerator (sum of fraction) of TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 4, 1, 3.2);
1283 h_cluster_phisize1_fraction_mean_numerator_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1284 h_cluster_phisize1_fraction_mean_numerator_side0[region]->GetYaxis()->SetTitle("Fraction");
1285 hm->registerHisto(h_cluster_phisize1_fraction_mean_numerator_side0[region]);
1286
1287 h_cluster_phisize1_fraction_mean_numerator_side1[region] = new TH1F(std::format("{}clusphisize1frac_mean_numerator_side1_{}", getHistoPrefix(), region).c_str(),
1288 std::format("Pt vs. Average fraction mean_numerator (sum of fraction) of TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 4, 1, 3.2);
1289 h_cluster_phisize1_fraction_mean_numerator_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1290 h_cluster_phisize1_fraction_mean_numerator_side1[region]->GetYaxis()->SetTitle("Fraction");
1291 hm->registerHisto(h_cluster_phisize1_fraction_mean_numerator_side1[region]);
1292
1293 h_cluster_phisize1_fraction_mean_denominator_side0[region] = new TH1F(std::format("{}clusphisize1frac_mean_denominator_side0_{}", getHistoPrefix(), region).c_str(),
1294 std::format("Pt vs. Average fraction mean_denominator (sum of track) of TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 4, 1, 3.2);
1295 h_cluster_phisize1_fraction_mean_denominator_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1296 h_cluster_phisize1_fraction_mean_denominator_side0[region]->GetYaxis()->SetTitle("Fraction");
1297 hm->registerHisto(h_cluster_phisize1_fraction_mean_denominator_side0[region]);
1298
1299 h_cluster_phisize1_fraction_mean_denominator_side1[region] = new TH1F(std::format("{}clusphisize1frac_mean_denominator_side1_{}", getHistoPrefix(), region).c_str(),
1300 std::format("Pt vs. Average fraction mean_denominator (sum of track) of TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 4, 1, 3.2);
1301 h_cluster_phisize1_fraction_mean_denominator_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1302 h_cluster_phisize1_fraction_mean_denominator_side1[region]->GetYaxis()->SetTitle("Fraction");
1303 hm->registerHisto(h_cluster_phisize1_fraction_mean_denominator_side1[region]);
1304 }
1305 }
1306
1307 std::pair<float, float> TpcSeedsQA::cal_tpc_eta_min_max(float vtxz)
1308 {
1309 float R = 780.;
1310 float HalfZ = 2110. / 2.;
1311 float theta_max = std::atan2(R, HalfZ - vtxz);
1312 float theta_min = std::atan2(R, -(vtxz + HalfZ));
1313 float eta_max = -std::log(std::tan(theta_max / 2));
1314 float eta_min = -std::log(std::tan(theta_min / 2));
1315 std::pair<float, float> min_max = std::make_pair(eta_min, eta_max);
1316 return min_max;
1317 }
1318
1319 float TpcSeedsQA::eta_to_theta(float eta)
1320 {
1321 return 2 * std::atan(std::exp(-eta));
1322 }
