File indexing completed on 2025-08-05 08:13:29
0001
0002
0003
0004
0005
0006
0007
0008 #include "AnaMvtxTestBeam2019.h"
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/PHTFileServer.h>
0012
0013 #include <phool/getClass.h>
0014
0015
0016
0017
0018
0019
0020 #include <trackbase/TrkrClusterv2.h>
0021 #include <trackbase/TrkrClusterContainerv1.h>
0022 #include <trackbase/TrkrHitSetContainerv1.h>
0023 #include <trackbase/TrkrHitSetv1.h>
0024
0025
0026 #include <TFile.h>
0027
0028 #include <TH1D.h>
0029 #include <TH2D.h>
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041 #define LogDebug(exp) std::cout<<"DEBUG: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0042 #define LogError(exp) std::cout<<"ERROR: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0043 #define LogWarning(exp) std::cout<<"WARNING: " <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0044
0045
0046
0047 AnaMvtxTestBeam2019::AnaMvtxTestBeam2019(const std::string &name,
0048 const std::string &ofName):
0049 SubsysReco(name),
0050 m_hits(nullptr),
0051 m_clusters(nullptr),
0052 m_lout_clusterQA(nullptr),
0053 m_lout_tracking(nullptr),
0054 m_foutname(ofName),
0055 do_tracking(false),
0056 m_ievent(0),
0057 m_ref_align_stave(1)
0058 {
0059 h1d_nevents = NULL;
0060 h1d_hit_layer = NULL;
0061 h1d_clu_layer = NULL;
0062 h1d_clu_map = NULL;
0063
0064 for (int il = 0; il < NLAYER; il++)
0065 {
0066 h2f_hit_map[il] = NULL;
0067 h2f_clu_xz[il] = NULL;
0068 h1d_clu_mean_dz[il] = NULL;
0069 h1d_clu_mean_dx[il] = NULL;
0070 h1d_clu_size[il] = NULL;
0071 h1d_clu_size_phi[il] = NULL;
0072 h1d_clu_size_z[il] = NULL;
0073
0074 htrk_dz[il] = NULL;
0075 htrk_dx[il] = NULL;
0076
0077 htrk_cut_dz[il] = NULL;
0078 htrk_cut_dx[il] = NULL;
0079 }
0080
0081 htrk = NULL;
0082 htrk_clus = NULL;
0083 htrk_chi2xy = NULL;
0084 htrk_chi2zy = NULL;
0085 htrk_cut = NULL;
0086 htrk_cut_clus = NULL;
0087 htrk_cut_chi2xy = NULL;
0088 htrk_cut_chi2zy = NULL;
0089
0090 m_mvtxtracking = new MvtxStandaloneTracking();
0091 }
0092
0093
0094
0095 int
0096 AnaMvtxTestBeam2019::Init(PHCompositeNode *topNode)
0097 {
0098 int chipColor[] = {kBlue, kRed, kGreen+2, kMagenta+2};
0099
0100
0101 m_lout_clusterQA = new TList();
0102 m_lout_clusterQA->SetName("ClusterQA");
0103 m_lout_clusterQA->SetOwner(true);
0104
0105
0106 h1d_nevents = new TH1D("h1d_nevents", ";;Nevts", 6, -.5, 5.5);
0107 h1d_hit_layer = new TH1D("h1d_hit_layer",
0108 "# hits/layer;layer;Counts", 6, -.5, 5.5);
0109 h1d_clu_layer = new TH1D("h1d_clu_layer",
0110 "# clusters/layer;layer;Counts", 6, -.5, 5.5);
0111 h1d_clu_map = new TH1D("h1d_clu_map",
0112 "cluster map;cluster_map;Counts", 16, -.5, 15.5);
0113 m_lout_clusterQA->Add(h1d_nevents);
0114 m_lout_clusterQA->Add(h1d_hit_layer);
0115 m_lout_clusterQA->Add(h1d_clu_layer);
0116 m_lout_clusterQA->Add(h1d_clu_map);
0117
0118 for (int il = 0; il < NLAYER; il++)
0119 {
0120 h2f_hit_map[il] = new TH2F(Form("h2d_hit_map_l%i", il),
0121 ";col [px];row [px]",
0122 9 * NCOL, -.5, (9 * NCOL) - .5,
0123 NROW, -.5, NROW - .5);
0124
0125 h2f_clu_xz[il] = new TH2F(Form("h2d_clu_xz_l%i", il),
0126 ";Z [cm];X [cm]",
0127 3E4, -15., 15.,
0128 2E3, -1., 1.);
0129
0130 h1d_clu_mean_dz[il] = new TH1D(Form("h1d_clu_mean_dz_l%i", il),
0131 ";#DeltaZ (cm);Counts/30 #mum",
0132 3E3, -1.5, 1.5);
0133 h1d_clu_mean_dz[il]->SetLineWidth(2);
0134 h1d_clu_mean_dz[il]->SetLineColor(chipColor[il]);
0135
0136 h1d_clu_mean_dx[il] = new TH1D(Form("h1d_clu_mean_dx_l%i", il),
0137 ";#DeltaX (cm);Counts/30 #mum",
0138 3E3, -1.5, 1.5);
0139 h1d_clu_mean_dx[il]->SetLineWidth(2);
0140 h1d_clu_mean_dx[il]->SetLineColor(chipColor[il]);
0141
0142 h1d_clu_size[il] = new TH1D(Form("h1d_clu_size_l%i", il),
0143 ";cluster size [px];Counts", 100, -.5, 99.5);
0144
0145 h1d_clu_size_phi[il] = new TH1D(Form("h1d_clu_size_phi_l%i", il),
0146 ";cluster size phi [px];Counts", 100, -.5, 99.5);
0147
0148 h1d_clu_size_z[il] = new TH1D(Form("h1d_clu_size_z_l%i", il),
0149 ";cluster size z [px];Counts", 100, -.5, 99.5);
0150
0151 m_lout_clusterQA->Add(h2f_hit_map[il]);
0152 m_lout_clusterQA->Add(h2f_clu_xz[il]);
0153 m_lout_clusterQA->Add(h1d_clu_mean_dz[il]);
0154 m_lout_clusterQA->Add(h1d_clu_mean_dx[il]);
0155 m_lout_clusterQA->Add(h1d_clu_size[il]);
0156 m_lout_clusterQA->Add(h1d_clu_size_phi[il]);
0157 m_lout_clusterQA->Add(h1d_clu_size_z[il]);
0158 }
0159
0160 if ( do_tracking )
0161 {
0162 m_lout_tracking = new TList();
0163 m_lout_tracking->SetName("Tracking");
0164 m_lout_tracking->SetOwner(true);
0165
0166
0167 htrk = new TH1D("htrk", ";trks / event", 16, -0.5, 15.5);
0168 htrk_clus = new TH1D("htrk_clus", ";#cluster/track;Counts", 10, -.5, 9.5);
0169 htrk_cut = new TH1D("htrk_cut", ";trks / event", 16, -0.5, 15.5);
0170 htrk_cut_clus = new TH1D("htrk_cut_clus", ";#cluster/track;Counts", 5, -.5, 4.5);
0171 m_lout_tracking->Add(htrk);
0172 m_lout_tracking->Add(htrk_clus);
0173 m_lout_tracking->Add(htrk_cut);
0174 m_lout_tracking->Add(htrk_cut_clus);
0175
0176 for (int il = 0; il < 4; il++)
0177 {
0178 htrk_dx[il] = new TH1D(Form("htrk_dx_l%i", il),
0179 ";track dx [cm]",
0180 2000, -.5, .5);
0181
0182 htrk_dz[il] = new TH1D(Form("htrk_dz_l%i", il),
0183 ";track dz [cm]",
0184 2000, -.5, .5);
0185
0186 htrk_cut_dx[il] = new TH1D(Form("htrk_cut_dx_l%i", il),
0187 ";track dx [pixels]",
0188 2000, -.5, .5);
0189
0190 htrk_cut_dz[il] = new TH1D(Form("htrk_cut_dz_l%i", il),
0191 ";track dz [cm]",
0192 2000, -.5, .5);
0193
0194 m_lout_tracking->Add(htrk_dx[il]);
0195 m_lout_tracking->Add(htrk_dz[il]);
0196 m_lout_tracking->Add(htrk_cut_dx[il]);
0197 m_lout_tracking->Add(htrk_cut_dz[il]);
0198 }
0199
0200 htrk_chi2xy = new TH1D("htrk_chi2xy",
0201 ";track chi2/ndf in x vs y",
0202 500, 0, 100);
0203
0204 htrk_chi2zy = new TH1D("htrk_chi2zy",
0205 ";track chi2/ndf in z vs y",
0206 500, 0, 100);
0207
0208 htrk_cut_chi2xy = new TH1D("htrk_cut_chi2xy",
0209 ";track chi2/ndf in x vs y",
0210 500, 0, 100);
0211
0212 htrk_cut_chi2zy = new TH1D("htrk_cut_chi2zy",
0213 ";track chi2/ndf in z vs y",
0214 500, 0, 100);
0215
0216 m_lout_tracking->Add(htrk_chi2xy);
0217 m_lout_tracking->Add(htrk_chi2zy);
0218 m_lout_tracking->Add(htrk_cut_chi2xy);
0219 m_lout_tracking->Add(htrk_cut_chi2zy);
0220 }
0221
0222 return Fun4AllReturnCodes::EVENT_OK;
0223 }
0224
0225
0226
0227 int
0228 AnaMvtxTestBeam2019::InitRun(PHCompositeNode *topNode)
0229 {
0230
0231 m_ievent = -1;
0232
0233 return Fun4AllReturnCodes::EVENT_OK;
0234 }
0235
0236
0237
0238 int
0239 AnaMvtxTestBeam2019::process_event(PHCompositeNode *topNode)
0240 {
0241 if ( Verbosity() > VERBOSITY_MORE )
0242 std::cout << "AnaMvtxTestBeam2019::process_event()" << std::endl;
0243
0244
0245 m_ievent++;
0246 h1d_nevents->Fill(0);
0247
0248
0249
0250
0251 int retval = GetNodes(topNode);
0252 if ( retval != Fun4AllReturnCodes::EVENT_OK )
0253 return retval;
0254
0255 h1d_nevents->Fill(1);
0256
0257
0258 unsigned int n_hits = 0;
0259 auto hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::mvtxId);
0260 for (auto hitsetitr = hitsetrange.first;
0261 hitsetitr != hitsetrange.second;
0262 ++hitsetitr)
0263 {
0264 auto hitset = hitsetitr->second;
0265
0266 auto layer = TrkrDefs::getLayer(hitset->getHitSetKey());
0267 auto stave = MvtxDefs::getStaveId(hitset->getHitSetKey());
0268 auto chip = MvtxDefs::getChipId(hitset->getHitSetKey());
0269 if (Verbosity() > VERBOSITY_A_LOT)
0270 {
0271 std::cout << "layer: " << layer << " stave: " \
0272 << stave << " chip: " << chip << std::endl;
0273 }
0274
0275
0276
0277
0278 auto hitrangei = hitset->getHits();
0279 for (auto hitr = hitrangei.first;
0280 hitr != hitrangei.second;
0281 ++hitr)
0282 {
0283 ++n_hits;
0284 auto col = MvtxDefs::getCol( hitr->first);
0285 auto row = MvtxDefs::getRow( hitr->first);
0286 auto stave_col = chip * SegmentationAlpide::NCols + col;
0287
0288 h1d_hit_layer->Fill(layer);
0289 h2f_hit_map[layer]->Fill(stave_col, row);
0290 }
0291 }
0292
0293 if (m_clusters->size() > n_hits)
0294 std::cout << __PRETTY_FUNCTION__ \
0295 << " : WARNING!! Nclus " << m_clusters->size() \
0296 << " > Nhits " << n_hits << std::endl;
0297
0298
0299
0300
0301 if (Verbosity() > VERBOSITY_MORE)
0302 {
0303 std::cout << "-- Looping over clusters." << std::endl;
0304 }
0305
0306 LyrClusMap m_LyrClusMap;
0307
0308 int n_clu_per_layer[NLAYER] = {0};
0309 float mean_clu_x[NLAYER] = {0.};
0310 float mean_clu_z[NLAYER] = {0.};
0311
0312 int n_clus = 0;
0313 auto clus_range = m_clusters->getClusters();
0314 for (auto iter_clus = clus_range.first;
0315 iter_clus != clus_range.second;
0316 ++iter_clus)
0317 {
0318 auto clus = iter_clus->second;
0319 auto ckey = clus->getClusKey();
0320 auto lyr = TrkrDefs::getLayer(ckey);
0321
0322 h1d_clu_layer->Fill(lyr);
0323
0324 if ((lyr < 0) || (lyr > 3)) continue;
0325
0326 h1d_clu_size[lyr]->Fill(clus->getAdc());
0327 h1d_clu_size_phi[lyr]->Fill(clus->getPhiSize()/SegmentationAlpide::PitchRow);
0328 h1d_clu_size_z[lyr]->Fill(clus->getZSize()/SegmentationAlpide::PitchCol);
0329
0330 h2f_clu_xz[lyr]->Fill(clus->getZ(), clus->getX());
0331
0332 n_clu_per_layer[lyr]++;
0333 mean_clu_z[lyr] += clus->getZ();
0334 mean_clu_x[lyr] += clus->getX();
0335
0336
0337 m_LyrClusMap.insert(std::make_pair(lyr, clus));
0338 n_clus++;
0339 }
0340
0341 char clu_map = 0;
0342 for (int ilyr = 0; ilyr < NLAYER; ++ilyr)
0343 {
0344 if (n_clu_per_layer[ilyr] > 0)
0345 {
0346 mean_clu_z[ilyr] /= (float)n_clu_per_layer[ilyr];
0347 mean_clu_x[ilyr] /= (float)n_clu_per_layer[ilyr];
0348
0349 h1d_clu_mean_dz[ilyr]->Fill(mean_clu_z[ilyr] - mean_clu_z[m_ref_align_stave]);
0350 h1d_clu_mean_dx[ilyr]->Fill(mean_clu_x[ilyr] - mean_clu_x[m_ref_align_stave]);
0351 clu_map |= ( 1 << ilyr );
0352 }
0353 }
0354 h1d_clu_map->Fill(clu_map);
0355
0356 if (Verbosity() > VERBOSITY_MORE)
0357 {
0358 std::cout << "-- evnt: " << m_ievent << std::endl;
0359 std::cout << "-- Nhits: " << n_hits << std::endl;
0360 std::cout << "-- Nclus: " << n_clus << std::endl;
0361 std::cout << "-- m_LyrClusMap.size(): " << m_LyrClusMap.size() << std::endl;
0362 }
0363
0364
0365
0366
0367 if (do_tracking)
0368 {
0369
0370
0371
0372 MvtxStandaloneTracking::MvtxTrackList tracklist;
0373
0374 std::vector<unsigned int> use_layers;
0375 use_layers.push_back(3);
0376 use_layers.push_back(2);
0377 use_layers.push_back(0);
0378 use_layers.push_back(1);
0379
0380 m_mvtxtracking->RunTracking(m_LyrClusMap, tracklist, use_layers);
0381
0382 htrk->Fill(tracklist.size());
0383
0384 int ncut = 0;
0385 for (unsigned int i = 0; i < tracklist.size(); i++)
0386 {
0387 htrk_clus->Fill(tracklist.at(i).ClusterList.size());
0388 htrk_chi2xy->Fill(tracklist.at(i).chi2_xy);
0389 htrk_chi2zy->Fill(tracklist.at(i).chi2_zy);
0390
0391 if (tracklist.at(i).chi2_xy < 5 && tracklist.at(i).chi2_zy < 5)
0392 {
0393 ++ncut;
0394 htrk_cut_chi2xy->Fill(tracklist.at(i).chi2_xy);
0395 htrk_cut_chi2zy->Fill(tracklist.at(i).chi2_zy);
0396 htrk_cut_clus->Fill(tracklist.at(i).ClusterList.size());
0397 }
0398
0399 for (unsigned int j = 0; j < tracklist.at(i).ClusterList.size(); j++)
0400 {
0401 auto ckey = tracklist.at(i).ClusterList.at(j)->getClusKey();
0402 auto lyr = TrkrDefs::getLayer(ckey);
0403
0404 if ((lyr < 0) || (lyr >= 4))
0405 {
0406 std::cout << PHWHERE << " WARNING: bad layer from track cluster. lyr:" << lyr << std::endl;
0407 continue;
0408 }
0409
0410 htrk_dx[lyr]->Fill(tracklist.at(i).dx.at(j));
0411 htrk_dz[lyr]->Fill(tracklist.at(i).dz.at(j));
0412
0413 if (tracklist.at(i).chi2_xy < 5 && tracklist.at(i).chi2_zy < 5)
0414 {
0415 htrk_cut_dx[lyr]->Fill(tracklist.at(i).dx.at(j));
0416 htrk_cut_dz[lyr]->Fill(tracklist.at(i).dz.at(j));
0417 }
0418 }
0419 }
0420 htrk_cut->Fill(ncut);
0421 }
0422
0423 return Fun4AllReturnCodes::EVENT_OK;
0424 }
0425
0426
0427
0428 int
0429 AnaMvtxTestBeam2019::End(PHCompositeNode *topNode)
0430 {
0431
0432 PHTFileServer::get().open(m_foutname, "RECREATE");
0433
0434 PHTFileServer::get().cd(m_foutname);
0435 gDirectory->mkdir(m_lout_clusterQA->GetName())->cd();
0436 m_lout_clusterQA->Write();
0437
0438 if (m_lout_tracking)
0439 {
0440 PHTFileServer::get().cd(m_foutname);
0441 gDirectory->mkdir(m_lout_tracking->GetName())->cd();
0442 m_lout_tracking->Write();
0443 }
0444 PHTFileServer::get().write(m_foutname);
0445 PHTFileServer::get().close();
0446
0447 return Fun4AllReturnCodes::EVENT_OK;
0448 }
0449
0450
0451
0452
0453
0454
0455 int
0456 AnaMvtxTestBeam2019::GetNodes(PHCompositeNode *topNode)
0457 {
0458
0459 m_hits = findNode::getClass<TrkrHitSetContainerv1>(topNode, "TRKR_HITSET");
0460 if (! m_hits)
0461 {
0462 std::cout << PHWHERE << " HITSETS node was not found on node tree" << std::endl;
0463 return Fun4AllReturnCodes::ABORTEVENT;
0464 }
0465
0466 m_clusters = findNode::getClass<TrkrClusterContainerv1>(topNode, "TRKR_CLUSTER");
0467 if (! m_clusters)
0468 {
0469 std::cout << PHWHERE << " TRKR_CLUSTER node not found on node tree" << std::endl;
0470 return Fun4AllReturnCodes::ABORTEVENT;
0471 }
0472
0473 return Fun4AllReturnCodes::EVENT_OK;
0474 }
0475
0476
0477
0478
0479
0480
0481 double
0482 AnaMvtxTestBeam2019::CalcSlope(double x0, double y0, double x1, double y1)
0483 {
0484 return (y1 - y0) / (x1 - x0);
0485 }
0486
0487
0488
0489
0490
0491
0492 double
0493 AnaMvtxTestBeam2019::CalcIntecept(double x0, double y0, double m)
0494 {
0495 return y0 - x0 * m;
0496 }
0497
0498
0499
0500
0501
0502
0503 double
0504 AnaMvtxTestBeam2019::CalcProjection(double x, double m, double b)
0505 {
0506 return m * x + b;
0507 }
