Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:21:27 

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/PHG4TpcGeom.h>
0015 #include <g4detectors/PHG4TpcGeomContainer.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<PHG4TpcGeomContainer>(topNode, m_g4GeomName);
0058   trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0059   vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vertexMapName);
0060 
0061   if (!trackmap || !clustermap || !actsgeom || !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 TPCGEOMCONTAINER" << std::endl;
0070     return Fun4AllReturnCodes::ABORTRUN;
0071   }
0072 
0073   // global position wrapper
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   // tracks with TPC clusters/tracklets
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   // h_trackcrossing_pos = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "trackcrossing_pos").c_str()));
0106   // h_trackcrossing_neg = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "trackcrossing_neg").c_str()));
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   // vertex
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   // h_vcrossing = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "vertexcrossing").c_str()));
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   // TPC has 3 regions, inner, mid and outer
0139   std::vector<int> region_layer_low = {7, 23, 39};
0140   std::vector<int> region_layer_high = {22, 38, 54};
0141 
0142   // make a layer to region multimap
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 &region : {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;  // the micromegas clusters are added to the TPC seeds
0188     }
0189     unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
0190     TrkrCluster *cluster = clustermap->findCluster(cluster_key);
0191     float adc = cluster->getAdc();
0192     PHG4TpcGeom *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     // Fully correct the cluster positions for the crossing and all distortions
0239     Acts::Vector3 global = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, track->get_crossing());
0240 
0241     // add the global positions to a vector to give to the cluster mover
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   // get the fully corrected cluster global positions
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     // Fully correct the cluster positions for the crossing and all distortions
0262     Acts::Vector3 global = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, track->get_crossing());
0263 
0264     // add the global positions to a vector to give to the cluster mover
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   // move the corrected cluster positions back to the original readout surface
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;  // the micromegas clusters are added to the TPC seeds
0316     }
0317 
0318     // only counts TPC R2 and R3
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;  // first: number of tracks not associated to a vertex, second: number of tracks associated to a vertex
0428   std::pair<int, int> ntrack_isfromvtx_neg;  // first: number of tracks not associated to a vertex, second: number of tracks associated to a vertex
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     // int trkcrossing = track->get_crossing();
0475 
0476     //    int nmaps = 0;
0477     //    int nintt = 0;
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       // case TrkrDefs::mvtxId:
0495       //   nmaps++;
0496       //   break;
0497       // case TrkrDefs::inttId:
0498       //   nintt++;
0499       //   break;
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       // std::cout << "No vertex found for track " << track->get_id() << std::std::endl;
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 
0561         if (std::fabs(eta) < 1.2)
0562         {
0563           h_ntpc_pos->Fill(ntpc);
0564         }
0565         h_ntpot_pos->Fill(nmms);
0566         h_ntpc_quality_pos->Fill(ntpc, quality);
0567         h_avgnclus_eta_phi_pos->Fill(eta, phi, ntpc);
0568         // h_trackcrossing_pos->Fill(trkcrossing);
0569         h_cluster_phisize1_fraction_pos->Fill((double) ntpc_phisize1 / (double) ntpc);
0570       }
0571     }
0572     else if (charge == -1)
0573     {
0574       h_ntpc_fullpt_neg->Fill(ntpc);
0575       if (dcapair_origin.second.first > 0)
0576       {
0577         h_dcaxyorigin_phi_north_neg->Fill(phi, dcapair_origin.first.first);
0578       }
0579       else if (dcapair_origin.second.first <= 0)
0580       {
0581         h_dcaxyorigin_phi_south_neg->Fill(phi, dcapair_origin.first.first);
0582       }
0583       h_dcazorigin_phi_neg->Fill(phi, dcapair_origin.second.first);
0584       if (pt > 1)
0585       {
0586         h_ntrack_neg->Fill(eta, phi);
0587         if (std::fabs(eta) < 1.2)
0588         {
0589           h_ntpc_neg->Fill(ntpc);
0590         }
0591         h_ntpot_neg->Fill(nmms);
0592         h_ntpc_quality_neg->Fill(ntpc, quality);
0593         h_avgnclus_eta_phi_neg->Fill(eta, phi, ntpc);
0594         // h_trackcrossing_neg->Fill(trkcrossing);
0595         h_cluster_phisize1_fraction_neg->Fill((double) ntpc_phisize1 / (double) ntpc);
0596       }
0597     }
0598   }
0599   h_ntrack1d_pos->Fill(ntrack1d_pos);
0600   h_ntrack1d_neg->Fill(ntrack1d_neg);
0601   h_ntrack1d_ptg1_pos->Fill(ntrack1d_ptg1_pos);
0602   h_ntrack1d_ptg1_neg->Fill(ntrack1d_ptg1_neg);
0603   h_ntrack1d_ptg1->Fill(ntrack1d_ptg1_pos + ntrack1d_ptg1_neg);
0604 
0605   h_ntrack_isfromvtx_pos->SetBinContent(1, h_ntrack_isfromvtx_pos->GetBinContent(1) + ntrack_isfromvtx_pos.first);
0606   h_ntrack_isfromvtx_pos->SetBinContent(2, h_ntrack_isfromvtx_pos->GetBinContent(2) + ntrack_isfromvtx_pos.second);
0607   h_ntrack_isfromvtx_neg->SetBinContent(1, h_ntrack_isfromvtx_neg->GetBinContent(1) + ntrack_isfromvtx_neg.first);
0608   h_ntrack_isfromvtx_neg->SetBinContent(2, h_ntrack_isfromvtx_neg->GetBinContent(2) + ntrack_isfromvtx_neg.second);
0609 
0610   // vertex
0611   h_nvertex->Fill(vertexmap->size());
0612   for (const auto &[key, vertex] : *vertexmap)
0613   {
0614     if (!vertex)
0615     {
0616       continue;
0617     }
0618 
0619     float vx = vertex->get_x();
0620     float vy = vertex->get_y();
0621     float vz = vertex->get_z();
0622     float vt = vertex->get_t0();
0623     float vchi2 = vertex->get_chisq();
0624     int vndof = vertex->get_ndof();
0625     // int vcrossing = vertex->get_beam_crossing();
0626 
0627     // std::cout << "vertex (x,y,z,t,chi2,ndof,crossing)=(" << vx << "," << vy << "," << vz << "," << vt << "," << vchi2 << "," << vndof << "," << vcrossing << ")" << std::endl;
0628 
0629     h_vx->Fill(vx);
0630     h_vy->Fill(vy);
0631     h_vx_vy->Fill(vx, vy);
0632     h_vz->Fill(vz);
0633     h_vt->Fill(vt);
0634     h_vchi2dof->Fill(float(vchi2 / vndof));
0635     // h_vcrossing->Fill(vcrossing);
0636 
0637     h_ntrackpervertex->Fill(vertex->size_tracks());
0638   }
0639 
0640   auto fill = [](TH1 *h, float val)
0641   { if (h) { h->Fill(val); } };
0642 
0643   std::set<unsigned int> tpc_seed_ids;
0644   for (const auto &[key, track] : *trackmap)
0645   {
0646     if (!track)
0647     {
0648       continue;
0649     }
0650     m_px = track->get_px();
0651     m_py = track->get_py();
0652     m_pz = track->get_pz();
0653     m_pt = std::sqrt((m_px * m_px) + (m_py * m_py));
0654     m_ptot = std::sqrt((m_px * m_px) + (m_py * m_py) + (m_pz * m_pz));
0655     TrackSeed *tpcseed = track->get_tpc_seed();
0656     m_charge = track->get_charge();
0657     m_dedx = calc_dedx(tpcseed);
0658 
0659     m_ntpc = 0;
0660     m_region.clear();
0661     m_clusgz.clear();
0662     m_cluslayer.clear();
0663     m_clusphisize.clear();
0664     m_cluszsize.clear();
0665 
0666     for (const auto &ckey : get_cluster_keys(track))
0667     {
0668       if (TrkrDefs::getTrkrId(ckey) == TrkrDefs::tpcId)
0669       {
0670         m_ntpc++;
0671       }
0672     }
0673 
0674     for (const auto &ckey : get_cluster_keys(track))
0675     {
0676       TrkrCluster *cluster = clustermap->findCluster(ckey);
0677       const Acts::Vector3 clusglob = m_globalPositionWrapper.getGlobalPositionDistortionCorrected(ckey, cluster, track->get_crossing());
0678 
0679       switch (TrkrDefs::getTrkrId(ckey))  // NOLINT(bugprone-switch-missing-default-case
0680       {
0681       case TrkrDefs::tpcId:
0682         const auto it = m_layerRegionMap.find(TrkrDefs::getLayer(ckey));
0683         int region = it->second;
0684         m_region.push_back(region);
0685         m_clusgz.push_back(clusglob.z());
0686         m_cluslayer.push_back(TrkrDefs::getLayer(ckey));
0687         m_clusphisize.push_back(cluster->getPhiSize());
0688         m_cluszsize.push_back(cluster->getZSize());
0689         int this_sector = (int) TpcDefs::getSectorId(ckey);
0690         int this_side = (int) TpcDefs::getSide(ckey);
0691         int is_onepad = 0;
0692         if (cluster->getPhiSize() <= 1)
0693         {
0694           is_onepad = 1;
0695         }
0696         if (m_ntpc > 30 && cluster->getPhiSize() > 1 && m_dedx < 1500 && m_pt > 1.0 && m_pt < 50)
0697         {
0698           h_adc_sector[region]->Fill((this_sector + 1) * (2 * (this_side - 0.5)), cluster->getAdc());
0699         }
0700         if (m_ntpc > 30 && m_dedx < 1000 && m_pt > 1.0 && m_pt < 50 && m_charge < 0)
0701         {
0702           h_onepad_frac[region]->Fill((this_sector + 1) * (2 * (this_side - 0.5)), is_onepad);
0703         }
0704         nclus[this_side][region][this_sector] += 1;
0705         madc[this_side][region][this_sector] += cluster->getAdc();
0706         break;
0707       }
0708     }
0709 
0710     if (m_ptot > 0.2 && m_ptot < 4 && m_ntpc > 30)
0711     {
0712       h_dedx->Fill(m_charge * m_ptot, m_dedx);
0713       if (collision_or_cosmics == false || (collision_or_cosmics == true && cal_track_length(track) > 25))
0714       {
0715         h_dedx_pcaz->Fill(track->get_z(), m_dedx);
0716       }
0717     }
0718 
0719     if (m_ptot > 1.0 && m_pt < 4 && m_ntpc > 30 && m_charge < 0 && m_dedx < 1000 && m_dedx > 50)
0720     {
0721       h_mip_dedx->Fill(m_dedx);
0722     }
0723 
0724     if (m_ntpc > 30)
0725     {
0726       float *cluster_dedx = cal_dedx_cluster(track);
0727       for (int iz = 0; iz < 10; iz++)
0728       {
0729         h_dedx_pq_z[iz]->Fill(m_charge * m_ptot, cluster_dedx[iz]);
0730       }
0731     }
0732 
0733     // if (m_pt > 1)
0734     //{
0735     //   h_ntpc->Fill(m_ntpc);
0736     // }
0737 
0738     int nClus = m_cluslayer.size();
0739     for (int cl = 0; cl < nClus; cl++)
0740     {
0741       if (m_pt > 1 && m_ntpc > 25)
0742       {
0743         if (m_clusphisize[cl] == 1 && m_cluszsize[cl] > 1)
0744         {
0745           if (m_clusgz[cl] < 0.)
0746           {
0747             const auto hiter = phihistos.find(m_region[cl]);
0748             if (hiter == phihistos.end())
0749             {
0750               continue;
0751             }
0752             fill(hiter->second.cphisize1pT_side0, m_pt);
0753           }
0754           else if (m_clusgz[cl] > 0.)
0755           {
0756             const auto hiter = phihistos.find(m_region[cl]);
0757             if (hiter == phihistos.end())
0758             {
0759               continue;
0760             }
0761             fill(hiter->second.cphisize1pT_side1, m_pt);
0762           }
0763         }
0764         if (m_clusphisize[cl] >= 1 && m_cluszsize[cl] > 1)
0765         {
0766           if (m_clusgz[cl] < 0.)
0767           {
0768             const auto hiter = phihistos.find(m_region[cl]);
0769             if (hiter == phihistos.end())
0770             {
0771               continue;
0772             }
0773             fill(hiter->second.cphisizegeq1pT_side0, m_pt);
0774           }
0775           else if (m_clusgz[cl] > 0.)
0776           {
0777             const auto hiter = phihistos.find(m_region[cl]);
0778             if (hiter == phihistos.end())
0779             {
0780               continue;
0781             }
0782             fill(hiter->second.cphisizegeq1pT_side1, m_pt);
0783           }
0784         }
0785       }
0786     }
0787 
0788     for (auto &pair : phihistos)
0789     {
0790       pair.second.Clear();
0791     }
0792 
0793     for (int cl = 0; cl < nClus; cl++)
0794     {
0795       if (m_pt > 1 && m_ntpc > 25)
0796       {
0797         if (m_clusgz[cl] < 0.)
0798         {
0799           const auto hiter = phihistos.find(m_region[cl]);
0800           if (hiter == phihistos.end())
0801           {
0802             continue;
0803           }
0804           if (m_clusphisize[cl] == 1 && m_cluszsize[cl] > 1)
0805           {
0806             hiter->second.ntpc_side0_phisize1++;
0807           }
0808           if (m_clusphisize[cl] >= 1 && m_cluszsize[cl] > 1)
0809           {
0810             hiter->second.ntpc_side0++;
0811           }
0812         }
0813         else if (m_clusgz[cl] > 0.)
0814         {
0815           const auto hiter = phihistos.find(m_region[cl]);
0816           if (hiter == phihistos.end())
0817           {
0818             continue;
0819           }
0820           if (m_clusphisize[cl] == 1 && m_cluszsize[cl] > 1)
0821           {
0822             hiter->second.ntpc_side1_phisize1++;
0823           }
0824           if (m_clusphisize[cl] >= 1 && m_cluszsize[cl] > 1)
0825           {
0826             hiter->second.ntpc_side1++;
0827           }
0828         }
0829       }
0830     }
0831 
0832     for (const auto &region : {0, 1, 2})
0833     {
0834       if (phihistos[region].ntpc_side0 > 0)
0835       {
0836         double frac_side0 = (double) phihistos[region].ntpc_side0_phisize1 / (double) phihistos[region].ntpc_side0;
0837         h_cluster_phisize1_fraction_side0[region]->Fill(frac_side0);
0838         h_cluster_phisize1_fraction_pt_side0[region]->Fill(m_pt, frac_side0);
0839 
0840         int index_pt_side0 = h_cluster_phisize1_fraction_mean_side0[region]->FindBin(m_pt) - 1;
0841         if (index_pt_side0 < h_cluster_phisize1_fraction_mean_side0[region]->GetNbinsX())
0842         {
0843           frac_side0_pt[region][index_pt_side0] += frac_side0;
0844           num_track_side0_pt[region][index_pt_side0]++;
0845         }
0846       }
0847 
0848       if (phihistos[region].ntpc_side1 > 0)
0849       {
0850         double frac_side1 = (double) phihistos[region].ntpc_side1_phisize1 / (double) phihistos[region].ntpc_side1;
0851         h_cluster_phisize1_fraction_side1[region]->Fill(frac_side1);
0852         h_cluster_phisize1_fraction_pt_side1[region]->Fill(m_pt, frac_side1);
0853 
0854         int index_pt_side1 = h_cluster_phisize1_fraction_mean_side1[region]->FindBin(m_pt) - 1;
0855         if (index_pt_side1 < h_cluster_phisize1_fraction_mean_side1[region]->GetNbinsX())
0856         {
0857           frac_side1_pt[region][index_pt_side1] += frac_side1;
0858           num_track_side1_pt[region][index_pt_side1]++;
0859         }
0860       }
0861     }
0862   }
0863 
0864   for (int iside = 0; iside < 2; iside++)
0865   {
0866     for (int iregion = 0; iregion < 3; iregion++)
0867     {
0868       for (int isector = 0; isector < 12; isector++)
0869       {
0870         //"event:segment:bco:side:region:sector:ncluster:meanadc"
0871         int nCluster = nclus[iside][iregion][isector];
0872         float meanAdc = 0;
0873         if (nCluster > 0)
0874         {
0875           meanAdc = madc[iside][iregion][isector] / (nCluster * 1.);
0876         }
0877         nt_sector_event_summary->Fill(m_event, m_segment, m_bco, iside, iregion, isector, nCluster, meanAdc);
0878       }
0879     }
0880   }
0881   m_event++;
0882   return Fun4AllReturnCodes::EVENT_OK;
0883 }
0884 
0885 std::vector<TrkrDefs::cluskey> TpcSeedsQA::get_cluster_keys(SvtxTrack *track)
0886 {
0887   std::vector<TrkrDefs::cluskey> out;
0888   for (const auto &seed : {track->get_silicon_seed(), track->get_tpc_seed()})
0889   {
0890     if (seed)
0891     {
0892       std::copy(seed->begin_cluster_keys(), seed->end_cluster_keys(), std::back_inserter(out));
0893     }
0894   }
0895   return out;
0896 }
0897 
0898 //____________________________________________________________________________..
0899 int TpcSeedsQA::EndRun(const int /*runnumber*/)
0900 {
0901   for (const auto &region : {0, 1, 2})
0902   {
0903     for (const auto &index_pt : {0, 1, 2, 3})
0904     {
0905       if (num_track_side0_pt[region][index_pt] > 0)
0906       {
0907         h_cluster_phisize1_fraction_mean_side0[region]->SetBinContent(index_pt + 1, frac_side0_pt[region][index_pt] / num_track_side0_pt[region][index_pt]);
0908         h_cluster_phisize1_fraction_mean_numerator_side0[region]->SetBinContent(index_pt + 1, frac_side0_pt[region][index_pt]);
0909         h_cluster_phisize1_fraction_mean_denominator_side0[region]->SetBinContent(index_pt + 1, num_track_side0_pt[region][index_pt]);
0910       }
0911       if (num_track_side1_pt[region][index_pt] > 0)
0912       {
0913         h_cluster_phisize1_fraction_mean_side1[region]->SetBinContent(index_pt + 1, frac_side1_pt[region][index_pt] / num_track_side1_pt[region][index_pt]);
0914         h_cluster_phisize1_fraction_mean_numerator_side1[region]->SetBinContent(index_pt + 1, frac_side1_pt[region][index_pt]);
0915         h_cluster_phisize1_fraction_mean_denominator_side1[region]->SetBinContent(index_pt + 1, num_track_side1_pt[region][index_pt]);
0916       }
0917     }
0918   }
0919 
0920   auto *hm = QAHistManagerDef::getHistoManager();
0921   assert(hm);
0922 
0923   return Fun4AllReturnCodes::EVENT_OK;
0924 }
0925 
0926 //____________________________________________________________________________..
0927 int TpcSeedsQA::End(PHCompositeNode * /*unused*/)
0928 {
0929   return Fun4AllReturnCodes::EVENT_OK;
0930 }
0931 
0932 std::string TpcSeedsQA::getHistoPrefix() const
0933 {
0934   return std::string("h_") + Name() + std::string("_");
0935 }
0936 
0937 void TpcSeedsQA::createHistos()
0938 {
0939   auto *hm = QAHistManagerDef::getHistoManager();
0940   assert(hm);
0941 
0942   {
0943     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);
0944     hm->registerHisto(h);
0945   }
0946 
0947   {
0948     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);
0949     hm->registerHisto(h);
0950   }
0951 
0952   {
0953     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);
0954     hm->registerHisto(h);
0955   }
0956 
0957   {
0958     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);
0959     hm->registerHisto(h);
0960   }
0961 
0962   {
0963     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);
0964     hm->registerHisto(h);
0965   }
0966 
0967   {
0968     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);
0969     hm->registerHisto(h);
0970   }
0971 
0972   {
0973     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);
0974     hm->registerHisto(h);
0975   }
0976 
0977   {
0978     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);
0979     hm->registerHisto(h);
0980   }
0981 
0982   {
0983     auto *h = new TH1F(std::string(getHistoPrefix() + "nrecotracks1d").c_str(), "Number of reconstructed tracks;Number of TPC tracklets;Entries", 50, 0, 200);
0984     hm->registerHisto(h);
0985   }
0986 
0987   {
0988     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);
0989     hm->registerHisto(h);
0990   }
0991 
0992   {
0993     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);
0994     hm->registerHisto(h);
0995   }
0996 
0997   {
0998     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);
0999     hm->registerHisto(h);
1000   }
1001 
1002   {
1003     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);
1004     hm->registerHisto(h);
1005   }
1006 
1007   {
1008     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);
1009     hm->registerHisto(h);
1010   }
1011 
1012   {
1013     auto *h = new TH1F(std::string(getHistoPrefix() + "pt").c_str(), "p_{T} distribution of reconstructed tracks;Track p_{T};Entries", 100, 0, 10);
1014     hm->registerHisto(h);
1015   }
1016 
1017   {
1018     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);
1019     hm->registerHisto(h);
1020   }
1021 
1022   {
1023     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);
1024     hm->registerHisto(h);
1025   }
1026 
1027   {
1028     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);
1029     hm->registerHisto(h);
1030   }
1031 
1032   {
1033     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);
1034     hm->registerHisto(h);
1035   }
1036 
1037   {
1038     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);
1039     hm->registerHisto(h);
1040   }
1041 
1042   {
1043     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);
1044     hm->registerHisto(h);
1045   }
1046 
1047   //  {
1048   //    auto h = new TH1F(std::string(getHistoPrefix() + "trackcrossing_pos").c_str(), "Positive track beam bunch crossing (pT>1GeV);Positive track crossing;Entries", 100, -100, 300);
1049   //    hm->registerHisto(h);
1050   //  }
1051 
1052   //  {
1053   //    auto h = new TH1F(std::string(getHistoPrefix() + "trackcrossing_neg").c_str(), "Negative track beam bunch crossing (pT>1GeV);Negative track crossing;Entries", 100, -100, 300);
1054   //    hm->registerHisto(h);
1055   //  }
1056 
1057   {
1058     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);
1059     hm->registerHisto(h);
1060   }
1061 
1062   {
1063     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);
1064     hm->registerHisto(h);
1065   }
1066 
1067   {
1068     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);
1069     hm->registerHisto(h);
1070   }
1071 
1072   {
1073     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);
1074     hm->registerHisto(h);
1075   }
1076 
1077   {
1078     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);
1079     hm->registerHisto(h);
1080   }
1081 
1082   {
1083     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);
1084     hm->registerHisto(h);
1085   }
1086 
1087   {
1088     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);
1089     hm->registerHisto(h);
1090   }
1091 
1092   {
1093     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);
1094     hm->registerHisto(h);
1095   }
1096 
1097   {
1098     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);
1099     hm->registerHisto(h);
1100   }
1101 
1102   {
1103     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);
1104     hm->registerHisto(h);
1105   }
1106 
1107   {
1108     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);
1109     hm->registerHisto(h);
1110   }
1111 
1112   {
1113     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);
1114     hm->registerHisto(h);
1115   }
1116 
1117   {
1118     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);
1119     hm->registerHisto(h);
1120   }
1121 
1122   {
1123     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);
1124     hm->registerHisto(h);
1125   }
1126 
1127   // vertex
1128   {
1129     auto *h = new TH1F(std::string(getHistoPrefix() + "nrecovertices").c_str(), "Num of reco vertices per event;Number of vertices;Entries", 20, 0, 20);
1130     hm->registerHisto(h);
1131   }
1132 
1133   {
1134     auto *h = new TH1F(std::string(getHistoPrefix() + "vx").c_str(), "Vertex x;Vertex x [cm];Entries", 100, -2.5, 2.5);
1135     hm->registerHisto(h);
1136   }
1137 
1138   {
1139     auto *h = new TH1F(std::string(getHistoPrefix() + "vy").c_str(), "Vertex y;Vertex y [cm];Entries", 100, -2.5, 2.5);
1140     hm->registerHisto(h);
1141   }
1142 
1143   {
1144     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);
1145     hm->registerHisto(h);
1146   }
1147 
1148   {
1149     auto *h = new TH1F(std::string(getHistoPrefix() + "vz").c_str(), "Vertex z;Vertex z [cm];Entries", 50, -25, 25);
1150     hm->registerHisto(h);
1151   }
1152 
1153   {
1154     auto *h = new TH1F(std::string(getHistoPrefix() + "vt").c_str(), "Vertex t;Vertex t [ns];Entries", 100, -1000, 20000);
1155     hm->registerHisto(h);
1156   }
1157 
1158   //  {
1159   //    auto h = new TH1F(std::string(getHistoPrefix() + "vertexcrossing").c_str(), "Vertex beam bunch crossing;Vertex crossing;Entries", 100, -100, 300);
1160   //    hm->registerHisto(h);
1161   //  }
1162 
1163   {
1164     auto *h = new TH1F(std::string(getHistoPrefix() + "vertexchi2dof").c_str(), "Vertex chi2/ndof;Vertex #chi2/ndof;Entries", 100, 0, 20);
1165     hm->registerHisto(h);
1166   }
1167 
1168   {
1169     auto *h = new TH1F(std::string(getHistoPrefix() + "ntrackspervertex").c_str(), "Num of tracks per vertex;Number of tracks per vertex;Entries", 50, 0, 50);
1170     hm->registerHisto(h);
1171   }
1172   {
1173     auto *h = new TH2F(std::string(getHistoPrefix() + "dedx").c_str(),
1174                        "Num of tracks per vertex;Number of tracks per vertex;Entries",
1175                        500, -2, 2, 500, 0, 3000);
1176     hm->registerHisto(h);
1177   }
1178 
1179   {
1180     auto *h = new TH1F(std::string(getHistoPrefix() + "mip_dedx").c_str(),
1181                        "dEdx of MIPs",
1182                        100, 0, 1000);
1183     hm->registerHisto(h);
1184   }
1185 
1186   {
1187     auto *h = new TH2F(std::string(getHistoPrefix() + "dedx_pcaz").c_str(),
1188                        "track dEdx vs pcaz",
1189                        100, -100, 100, 50, 0, 2000);
1190     hm->registerHisto(h);
1191   }
1192 
1193   for (int i = 0; i < 10; i++)
1194   {
1195     h_dedx_pq_z[i] = new TH2F(std::format("{}dedx_pq_{}", getHistoPrefix(), i).c_str(),
1196                               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);
1197     hm->registerHisto(h_dedx_pq_z[i]);
1198   }
1199 
1200   {
1201     auto *nt = new TNtuple(std::string(getHistoPrefix() + "sector_event_summary").c_str(),
1202                            "sector_event_summary", "event:segment:bco:side:region:sector:ncluster:meanadc");
1203     hm->registerHisto(nt);
1204   }
1205 
1206   for (const auto &region : {0, 1, 2})
1207   {
1208     h_adc_sector[region] = new TH2F(std::format("{}adc_sector_{}", getHistoPrefix(), region).c_str(),
1209                                     std::format("ADC spectrum per, region_{}", region).c_str(), 25, -12.5, 12.5, 50, 0, 1500);
1210     hm->registerHisto(h_adc_sector[region]);
1211 
1212     h_onepad_frac[region] = new TProfile(std::format("{}onepad_frac_{}", getHistoPrefix(), region).c_str(),
1213                                          std::format("TPC Cluster Phi Size == 1 fraction per sector, region_{}", region).c_str(), 25, -12.5, 12.5);
1214     hm->registerHisto(h_onepad_frac[region]);
1215 
1216     h_clusphisize1pt_side0[region] = new TH1F(std::format("{}clusphisize1pT_side0_{}", getHistoPrefix(), region).c_str(),
1217                                               std::format("TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 4, 1, 3.2);
1218     h_clusphisize1pt_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1219     hm->registerHisto(h_clusphisize1pt_side0[region]);
1220 
1221     h_clusphisize1pt_side1[region] = new TH1F(std::format("{}clusphisize1pT_side1_{}", getHistoPrefix(), region).c_str(),
1222                                               std::format("TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 4, 1, 3.2);
1223     h_clusphisize1pt_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1224     hm->registerHisto(h_clusphisize1pt_side1[region]);
1225 
1226     h_clusphisizegeq1pt_side0[region] = new TH1F(std::format("{}clusphisizegeq1pT_side0_{}", getHistoPrefix(), region).c_str(),
1227                                                  std::format("TPC Cluster Phi Size >= 1, side 0, region_{}", region).c_str(), 4, 1, 3.2);
1228     h_clusphisizegeq1pt_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1229     hm->registerHisto(h_clusphisizegeq1pt_side0[region]);
1230 
1231     h_clusphisizegeq1pt_side1[region] = new TH1F(std::format("{}clusphisizegeq1pT_side1_{}", getHistoPrefix(), region).c_str(),
1232                                                  std::format("TPC Cluster Phi Size >= 1, side 1, region_{}", region).c_str(), 4, 1, 3.2);
1233     h_clusphisizegeq1pt_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1234     hm->registerHisto(h_clusphisizegeq1pt_side1[region]);
1235 
1236     h_cluster_phisize1_fraction_side0[region] = new TH1F(std::format("{}clusphisize1frac_side0_{}", getHistoPrefix(), region).c_str(),
1237                                                          std::format("Fraction of TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 100, 0, 1);
1238     h_cluster_phisize1_fraction_side0[region]->GetXaxis()->SetTitle("Fraction");
1239     hm->registerHisto(h_cluster_phisize1_fraction_side0[region]);
1240 
1241     h_cluster_phisize1_fraction_side1[region] = new TH1F(std::format("{}clusphisize1frac_side1_{}", getHistoPrefix(), region).c_str(),
1242                                                          std::format("Fraction of TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 100, 0, 1);
1243     h_cluster_phisize1_fraction_side1[region]->GetXaxis()->SetTitle("Fraction");
1244     hm->registerHisto(h_cluster_phisize1_fraction_side1[region]);
1245 
1246     h_cluster_phisize1_fraction_pt_side0[region] = new TH2F(std::format("{}clusphisize1frac_pt_side0_{}", getHistoPrefix(), region).c_str(),
1247                                                             std::format("Pt vs. Fraction of TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 4, 1, 3.2, 100, 0, 1);
1248     h_cluster_phisize1_fraction_pt_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1249     h_cluster_phisize1_fraction_pt_side0[region]->GetYaxis()->SetTitle("Fraction");
1250     hm->registerHisto(h_cluster_phisize1_fraction_pt_side0[region]);
1251 
1252     h_cluster_phisize1_fraction_pt_side1[region] = new TH2F(std::format("{}clusphisize1frac_pt_side1_{}", getHistoPrefix(), region).c_str(),
1253                                                             std::format("Pt vs. Fraction of TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 4, 1, 3.2, 100, 0, 1);
1254     h_cluster_phisize1_fraction_pt_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1255     h_cluster_phisize1_fraction_pt_side1[region]->GetYaxis()->SetTitle("Fraction");
1256     hm->registerHisto(h_cluster_phisize1_fraction_pt_side1[region]);
1257 
1258     h_cluster_phisize1_fraction_mean_side0[region] = new TH1F(std::format("{}clusphisize1frac_mean_side0_{}", getHistoPrefix(), region).c_str(),
1259                                                               std::format("Pt vs. Average fraction of TPC Cluster Phi Size == 1, side 0, region_{}", region).c_str(), 4, 1, 3.2);
1260     h_cluster_phisize1_fraction_mean_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1261     h_cluster_phisize1_fraction_mean_side0[region]->GetYaxis()->SetTitle("Fraction");
1262     hm->registerHisto(h_cluster_phisize1_fraction_mean_side0[region]);
1263 
1264     h_cluster_phisize1_fraction_mean_side1[region] = new TH1F(std::format("{}clusphisize1frac_mean_side1_{}", getHistoPrefix(), region).c_str(),
1265                                                               std::format("Pt vs. Average fraction of TPC Cluster Phi Size == 1, side 1, region_{}", region).c_str(), 4, 1, 3.2);
1266     h_cluster_phisize1_fraction_mean_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1267     h_cluster_phisize1_fraction_mean_side1[region]->GetYaxis()->SetTitle("Fraction");
1268     hm->registerHisto(h_cluster_phisize1_fraction_mean_side1[region]);
1269 
1270     h_cluster_phisize1_fraction_mean_numerator_side0[region] = new TH1F(std::format("{}clusphisize1frac_mean_numerator_side0_{}", getHistoPrefix(), region).c_str(),
1271                                                                         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);
1272     h_cluster_phisize1_fraction_mean_numerator_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1273     h_cluster_phisize1_fraction_mean_numerator_side0[region]->GetYaxis()->SetTitle("Fraction");
1274     hm->registerHisto(h_cluster_phisize1_fraction_mean_numerator_side0[region]);
1275 
1276     h_cluster_phisize1_fraction_mean_numerator_side1[region] = new TH1F(std::format("{}clusphisize1frac_mean_numerator_side1_{}", getHistoPrefix(), region).c_str(),
1277                                                                         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);
1278     h_cluster_phisize1_fraction_mean_numerator_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1279     h_cluster_phisize1_fraction_mean_numerator_side1[region]->GetYaxis()->SetTitle("Fraction");
1280     hm->registerHisto(h_cluster_phisize1_fraction_mean_numerator_side1[region]);
1281 
1282     h_cluster_phisize1_fraction_mean_denominator_side0[region] = new TH1F(std::format("{}clusphisize1frac_mean_denominator_side0_{}", getHistoPrefix(), region).c_str(),
1283                                                                           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);
1284     h_cluster_phisize1_fraction_mean_denominator_side0[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1285     h_cluster_phisize1_fraction_mean_denominator_side0[region]->GetYaxis()->SetTitle("Fraction");
1286     hm->registerHisto(h_cluster_phisize1_fraction_mean_denominator_side0[region]);
1287 
1288     h_cluster_phisize1_fraction_mean_denominator_side1[region] = new TH1F(std::format("{}clusphisize1frac_mean_denominator_side1_{}", getHistoPrefix(), region).c_str(),
1289                                                                           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);
1290     h_cluster_phisize1_fraction_mean_denominator_side1[region]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
1291     h_cluster_phisize1_fraction_mean_denominator_side1[region]->GetYaxis()->SetTitle("Fraction");
1292     hm->registerHisto(h_cluster_phisize1_fraction_mean_denominator_side1[region]);
1293   }
1294 }
1295 
1296 std::pair<float, float> TpcSeedsQA::cal_tpc_eta_min_max(float vtxz)
1297 {
1298   float R = 780.;
1299   float HalfZ = 2110. / 2.;
1300   float theta_max = std::atan2(R, HalfZ - vtxz);
1301   float theta_min = std::atan2(R, -(vtxz + HalfZ));
1302   float eta_max = -std::log(std::tan(theta_max / 2));
1303   float eta_min = -std::log(std::tan(theta_min / 2));
1304   std::pair<float, float> min_max = std::make_pair(eta_min, eta_max);
1305   return min_max;
1306 }
1307 
1308 float TpcSeedsQA::eta_to_theta(float eta)
1309 {
1310   return 2 * std::atan(std::exp(-eta));
1311 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!