QA/MvtxMatchingEfficiencyWithShapes/MvtxMatchingEfficiencyWithShapes.cc

0001 #include "MvtxMatchingEfficiencyWithShapes.h"
0002
0003 #include <g4detectors/PHG4TpcCylinderGeom.h>
0004 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0005
0006 #include <trackbase/ActsGeometry.h>
0007 #include <trackbase/InttDefs.h>
0008 #include <trackbase/MvtxDefs.h>
0009 #include <trackbase/TrkrCluster.h>
0010 #include <trackbase/TrkrClusterContainer.h>
0011 #include <trackbase/TrkrClusterHitAssoc.h>
0012 #include <trackbase/TrkrDefs.h>
0013 #include <trackbase/TrkrHitSetContainer.h>
0014
0015 #include <trackbase_historic/SvtxTrackMap.h>
0016 #include <trackbase_historic/TrackAnalysisUtils.h>
0017 #include <trackbase_historic/TrackSeed.h>
0018 #include <trackbase_historic/TrackSeedContainer.h>
0019 #include <trackbase_historic/SvtxTrack.h>              // for SvtxTrack
0020 #include <trackbase_historic/SvtxTrackState.h>         // for SvtxTrackState
0021
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023 #include <fun4all/PHTFileServer.h>
0024 #include <fun4all/SubsysReco.h>
0025
0026 #include <phool/PHCompositeNode.h>
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h>
0029
0030 #include <TH1.h>
0031 #include <TTree.h>
0032
0033 #include <algorithm>
0034 #include <cmath>
0035 #include <iostream>                                    // for basic_ostream
0036 #include <set>
0037 #include <utility>                                     // for get, pair
0038
0039 namespace
0040 {
0041   //! range adaptor to be able to use range-based for loop
0042   template <class T>
0043   class range_adaptor
0044   {
0045    public:
0046     explicit range_adaptor(const T& range)
0047       : m_range(range)
0048     {
0049     }
0050     const typename T::first_type& begin() { return m_range.first; }
0051     const typename T::second_type& end() { return m_range.second; }
0052
0053    private:
0054     T m_range;
0055   };
0056 }  // namespace
0057
0058 //____________________________________________________________________________..
0059 MvtxMatchingEfficiencyWithShapes::MvtxMatchingEfficiencyWithShapes(const std::string& name, const std::string& outputfilename)
0060   : SubsysReco(name)
0061   , m_outputFileName(outputfilename)
0062 {
0063 }
0064
0065 //____________________________________________________________________________..
0066 int MvtxMatchingEfficiencyWithShapes::InitRun(PHCompositeNode* /*topNode*/)
0067 {
0068   // print configuration
0069   std::cout << "MvtxMatchingEfficiencyWithShapes::InitRun " << std::endl;
0070
0071   PHTFileServer::open(m_outputFileName, "RECREATE");
0072
0073   h_status = new TH1I("h_status_wTPC", "h_status_wTPC", 8, 0, 8);
0074   h_status->GetXaxis()->SetBinLabel(1, "Total");
0075   h_status->GetXaxis()->SetBinLabel(2, "INTT >= 2");
0076   h_status->GetXaxis()->SetBinLabel(3, "INTT same crossing");
0077   h_status->GetXaxis()->SetBinLabel(4, "INTT unique clusters");
0078   h_status->GetXaxis()->SetBinLabel(5, "MVTX unique clusters");
0079   h_status->GetXaxis()->SetBinLabel(6, "z < 10 cm");
0080   h_status->GetXaxis()->SetBinLabel(7, "eta < 1.1");
0081   h_status->GetXaxis()->SetBinLabel(8, "MVTX >= 2 layer");
0082
0083   h_INTT_time_delta = new TH1I("h_INTT_time_delta_wTPC", "h_INTT_time_delta_wTPC", 900, 0.5, 900.5);
0084
0085   // Create a new TTree
0086   tree = new TTree("TPC_tracks", "TPC_tracks");
0087
0088   // Create branches
0089   tree->Branch("pt", &pt, "pt/F");
0090   tree->Branch("eta", &eta, "eta/F");
0091   tree->Branch("phi", &phi, "phi/F");
0092   tree->Branch("frac_p_z", &frac_p_z, "frac_p_z/F");
0093   tree->Branch("dEdx", &dEdx, "dEdx/F");
0094   tree->Branch("nTPC", &nTPC, "nTPC/I");
0095   tree->Branch("layers", &layers, "layers/I");
0096   tree->Branch("states", &states, "states/I");
0097 // NOLINTNEXTLINE (readability-container-data-pointer)
0098   tree->Branch("shape_L0_C0_nhits", nhits_arr[0].data(), "shape_L0_C0_nhits/I");
0099   tree->Branch("shape_L0_C0_x", x_arr[0].data(), "shape_L0_C0_x/I");
0100   tree->Branch("shape_L0_C0_y", y_arr[0].data(), "shape_L0_C0_y/I");
0101   tree->Branch("shape_L0_C0_key", key_arr[0].data(), "shape_L0_C0_key/l");
0102
0103   tree->Branch("shape_L0_C1_nhits", &nhits_arr[0][1], "shape_L0_C1_nhits/I");
0104   tree->Branch("shape_L0_C1_x", &x_arr[0][1], "shape_L0_C1_x/I");
0105   tree->Branch("shape_L0_C1_y", &y_arr[0][1], "shape_L0_C1_y/I");
0106   tree->Branch("shape_L0_C1_key", &key_arr[0][1], "shape_L0_C1_key/l");
0107
0108   tree->Branch("shape_L1_C0_nhits", nhits_arr[1].data(), "shape_L1_C0_nhits/I");
0109   tree->Branch("shape_L1_C0_x", x_arr[1].data(), "shape_L1_C0_x/I");
0110   tree->Branch("shape_L1_C0_y", y_arr[1].data(), "shape_L1_C0_y/I");
0111   tree->Branch("shape_L1_C0_key", key_arr[1].data(), "shape_L1_C0_key/l");
0112
0113   tree->Branch("shape_L1_C1_nhits", &nhits_arr[1][1], "shape_L1_C1_nhits/I");
0114   tree->Branch("shape_L1_C1_x", &x_arr[1][1], "shape_L1_C1_x/I");
0115   tree->Branch("shape_L1_C1_y", &y_arr[1][1], "shape_L1_C1_y/I");
0116   tree->Branch("shape_L1_C1_key", &key_arr[1][1], "shape_L1_C1_key/l");
0117
0118   tree->Branch("shape_L2_C0_nhits", nhits_arr[2].data(), "shape_L2_C0_nhits/I");
0119   tree->Branch("shape_L2_C0_x", x_arr[2].data(), "shape_L2_C0_x/I");
0120   tree->Branch("shape_L2_C0_y", y_arr[2].data(), "shape_L2_C0_y/I");
0121   tree->Branch("shape_L2_C0_key", key_arr[2].data(), "shape_L2_C0_key/l");
0122
0123   tree->Branch("shape_L2_C2_nhits", &nhits_arr[2][1], "shape_L2_C2_nhits/I");
0124   tree->Branch("shape_L2_C2_x", &x_arr[2][1], "shape_L2_C2_x/I");
0125   tree->Branch("shape_L2_C2_y", &y_arr[2][1], "shape_L2_C2_y/I");
0126   tree->Branch("shape_L2_C2_key", &key_arr[2][1], "shape_L2_C2_key/l");
0127
0128   return Fun4AllReturnCodes::EVENT_OK;
0129 }
0130
0131 //____________________________________________________________________________..
0132 int MvtxMatchingEfficiencyWithShapes::process_event(PHCompositeNode* topNode)
0133 {
0134   ievent++;
0135
0136   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0137   if (!trackmap)
0138   {
0139     std::cout << PHWHERE
0140               << "SvtxTrackMap node is missing, can't collect particles"
0141               << std::endl;
0142     return -1;
0143   }
0144
0145   TrackSeedContainer* silicon_track_map = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0146   if (!silicon_track_map)
0147   {
0148     std::cout << PHWHERE
0149               << "SiliconTrackSeedContainer node is missing, can't collect particles"
0150               << std::endl;
0151     return -1;
0152   }
0153
0154   m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0155   if (!m_hitsetcontainer)
0156   {
0157     std::cout << "m_hitsetcontainer not found" << std::endl;
0158   }
0159   // assert(m_hitsetcontainer);
0160
0161   // cluster hit association map
0162   m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0163   if (!m_cluster_hit_map)
0164   {
0165     std::cout << "m_cluster_hit_map not found" << std::endl;
0166   }
0167
0168   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0169   if (!_cluster_map)
0170   {
0171     _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0172   }
0173
0174   if (!_cluster_map)
0175   {
0176     std::cout << PHWHERE
0177               << "TrkrClusterContainer node is missing"
0178               << std::endl;
0179     return -1;
0180   }
0181
0182   _geom_container = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0183   if (!_geom_container)
0184   {
0185     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0186     return Fun4AllReturnCodes::ABORTEVENT;
0187   }
0188
0189   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0190   if (!m_tGeometry)
0191   {
0192     std::cout << PHWHERE << "No acts reco geometry, bailing.";
0193     return Fun4AllReturnCodes::ABORTEVENT;
0194   }
0195
0196   // std::cout<<"SvtxTrackMap"<<std::endl;
0197   // fetch MVTX clusters
0198   std::vector<TrkrDefs::cluskey> cluster_vector;
0199   std::vector<TrkrDefs::cluskey> cluster_vector_INTT;
0200   for (auto& it : *trackmap)
0201   {
0202     for (const auto& ckey : TrackAnalysisUtils::get_cluster_keys(it.second))
0203     {
0204       switch (TrkrDefs::getTrkrId(ckey))
0205       {
0206       case TrkrDefs::mvtxId:
0207         cluster_vector.push_back(ckey);
0208         break;
0209       case TrkrDefs::inttId:
0210         cluster_vector_INTT.push_back(ckey);
0211         break;
0212       default:
0213     break;
0214       }
0215     }
0216   }
0217
0218   // find duplicated MVTX clusters in an event
0219   std::set<TrkrDefs::cluskey> duplicate_cluster_vector = findDuplicates(cluster_vector);
0220   // find duplicated INTT clusters in an event
0221   std::set<TrkrDefs::cluskey> duplicate_cluster_vector_INTT = findDuplicates(cluster_vector_INTT);
0222
0223   for (auto& it : *trackmap)
0224   {
0225     // number of cluster per subsystem
0226     int n_intt_hits = 0;
0227     int n_tpc_hits = 0;
0228     // fired mvtx layers 0 - innermost , 2 - outermost
0229     bool mvtx_l[3] = {false, false, false};
0230     bool mvtx_l_state[3] = {false, false, false};
0231     std::vector<float> intt_time;
0232
0233     SvtxTrack* track = it.second;
0234     // count clusters per subsystem and fired mvtx layers
0235     for (const auto& ckey : TrackAnalysisUtils::get_cluster_keys(track))
0236     {
0237       switch (TrkrDefs::getTrkrId(ckey))
0238       {
0239       case TrkrDefs::mvtxId:
0240         mvtx_l[TrkrDefs::getLayer(ckey)] = true;
0241         break;
0242       case TrkrDefs::inttId:
0243         intt_time.push_back(InttDefs::getTimeBucketId(ckey));
0244         n_intt_hits++;
0245         break;
0246       case TrkrDefs::tpcId:
0247         n_tpc_hits++;
0248         break;
0249       default:
0250     break;
0251       }
0252     }
0253
0254     // check states in MVTX layers
0255     for (auto state_it = track->begin_states(); state_it != track->end_states(); ++state_it)
0256     {
0257       auto clus_key = state_it->second->get_cluskey();
0258       switch (TrkrDefs::getTrkrId(clus_key))
0259       {
0260       case TrkrDefs::mvtxId:
0261         mvtx_l_state[TrkrDefs::getLayer(clus_key)] = true;
0262         break;
0263       default:
0264     break;
0265       }
0266     }
0267
0268     // all events
0269     h_status->Fill(0.5);
0270
0271     // require number of INTT clusters >= 2
0272     if (n_intt_hits < 2)
0273     {
0274       continue;
0275     }
0276     h_status->Fill(1.5);
0277
0278     // require the INTT clusters are from same crossing
0279     bool good_intt_time = std::all_of(intt_time.begin() + 1, intt_time.end(), [&](int i)
0280                                       { return i == intt_time[0]; });
0281
0282     if (good_intt_time == false)
0283     {
0284       if (n_intt_hits == 2)
0285       {
0286         h_INTT_time_delta->Fill(std::abs(intt_time[0] - intt_time[1]));
0287       }
0288       continue;
0289     }
0290     intt_time.clear();
0291     h_status->Fill(2.5);
0292
0293     // require the tracks does not share INTT clusters
0294     bool good_track_INTT = true;
0295     for (const auto& ckey : TrackAnalysisUtils::get_cluster_keys(track))
0296     {
0297       if (duplicate_cluster_vector_INTT.contains(ckey))
0298       {
0299         good_track_INTT = false;
0300       }
0301     }
0302     if (good_track_INTT == false)
0303     {
0304       continue;
0305     }
0306     h_status->Fill(3.5);
0307
0308     // require the tracks does not share MVTX clusters
0309     bool good_track = true;
0310     for (const auto& ckey : TrackAnalysisUtils::get_cluster_keys(track))
0311     {
0312       if (duplicate_cluster_vector.contains(ckey))
0313       {
0314         good_track = false;
0315       }
0316     }
0317     if (good_track == false)
0318     {
0319       continue;
0320     }
0321     h_status->Fill(4.5);
0322
0323     // z position cut
0324     if (std::abs(track->get_z()) > 10)
0325     {
0326       continue;
0327     }
0328     h_status->Fill(5.5);
0329
0330     // eta cut
0331     if (std::abs(track->get_eta()) > 1.1)
0332     {
0333       continue;
0334     }
0335     h_status->Fill(6.5);
0336
0337     // require 2 or 3 MVTX layers fired (layers not clusters)
0338     if (std::count_if(std::begin(mvtx_l), std::end(mvtx_l), [](bool i)
0339                       { return i == true; }) < 2)
0340     {
0341       continue;
0342     }
0343     h_status->Fill(7.5);
0344
0345     // fill Ttree variables
0346     pt = track->get_pt();
0347     eta = track->get_eta();
0348     phi = track->get_phi();
0349     frac_p_z = track->get_p() / track->get_charge();
0350     dEdx = calc_dedx(track->get_tpc_seed());
0351     layers = -1;
0352     states = -1;
0353     nTPC = n_tpc_hits;
0354
0355     // initiate cluster shape variables
0356     int layer_cluster_counter[3] = {-1, -1, -1};
0357     for (auto& inner_array : nhits_arr)
0358     {
0359       std::fill(inner_array.begin(), inner_array.end(), 0);
0360     }
0361     for (auto& inner_array : x_arr)
0362     {
0363       std::fill(inner_array.begin(), inner_array.end(), 0);
0364     }
0365     for (auto& inner_array : y_arr)
0366     {
0367       std::fill(inner_array.begin(), inner_array.end(), 0);
0368     }
0369     for (auto& inner_array : key_arr)
0370     {
0371       std::fill(inner_array.begin(), inner_array.end(), 0);
0372     }
0373
0374     // calculate cluster shape
0375     for (const auto& ckey : TrackAnalysisUtils::get_cluster_keys(track))
0376     {
0377       switch (TrkrDefs::getTrkrId(ckey))
0378       {
0379       case TrkrDefs::mvtxId:
0380       {
0381         int layer = TrkrDefs::getLayer(ckey);
0382         // std::cout << "layer " << layer << std::endl;
0383         layer_cluster_counter[layer]++;
0384         int cluster_index = layer_cluster_counter[layer];
0385
0386         // find associated hits
0387         const auto hit_range = m_cluster_hit_map->getHits(ckey);
0388
0389         int cL = 2000;
0390         int cR = 0;
0391         int cT = 2000;
0392         int cB = 0;
0393         uint64_t cmap = 0;
0394
0395         for (const auto& [ckey2, hitkey] : range_adaptor(hit_range))
0396         {
0397           int col = MvtxDefs::getCol(hitkey);
0398           int row = MvtxDefs::getRow(hitkey);
0399           cL = std::min(cL, col);
0400           cR = std::max(cR, col);
0401           cT = std::min(cT, row);
0402           cB = std::max(cB, row);
0403         }
0404
0405         int x = cR - cL + 1;
0406         int y = cB - cT + 1;
0407         int nhits = std::distance(hit_range.first, hit_range.second);
0408
0409         // key is saved only for shapes < 64 hits
0410         if (x * y < 64)
0411         {
0412           for (const auto& [ckey2, hitkey] : range_adaptor(hit_range))
0413           {
0414         // NOLINTNEXTLINE(hicpp-signed-bitwise)
0415             cmap |= (1ULL << ((MvtxDefs::getCol(hitkey) - cL) + ((cR - cL + 1) * (MvtxDefs::getRow(hitkey) - cT))));
0416             // std::cout<<"toggle bit "<<((MvtxDefs::getCol(hitkey)-cL)+((cR-cL+1)*(MvtxDefs::getRow(hitkey)-cT)))<<" "<< (MvtxDefs::getCol(hitkey)-cL)<<" "<<(cR-cL+1) <<" "<< (-MvtxDefs::getRow(hitkey)+cT)<<std::endl;
0417           }
0418           // std::cout<<"key "<<cmap<<std::endl;
0419         }
0420
0421         // shape debug prints
0422
0423         // std::cout<<"L: "<<cL<<" R: "<<cR<<" B: "<<cB<<" T: "<<cT<<std::endl;
0424         // std::cout<<"key int: "<<cmap<<std::endl;
0425         // std::cout<<"key uint: "<<(unsigned long long)cmap<<std::endl;
0426         // std::cout<<"key int: ";
0427
0428         // std::bitset<64> bits(cmap);
0429         // std::cout << bits << std::endl;
0430
0431         // Using bitwise operations
0432         // for (int i = 63; i >= 0; --i) {
0433         //    std::cout << ((cmap >> i) & 1);
0434         //}
0435         // std::cout << std::endl;
0436
0437         // for( const auto& [ckey2, hitkey]:range_adaptor( hit_range ) ){
0438         //   std::cout<<"row: "<<MvtxDefs::getRow(hitkey)<<" col "<< MvtxDefs::getCol(hitkey)<<" bitmap "<<((MvtxDefs::getCol(hitkey)-cL)+((cR-cL+1)*(MvtxDefs::getRow(hitkey)-cT)))<<std::endl;
0439         // }
0440
0441         // Fill the arrays
0442         nhits_arr[layer][cluster_index] = nhits;
0443         x_arr[layer][cluster_index] = x;
0444         y_arr[layer][cluster_index] = y;
0445         key_arr[layer][cluster_index] = cmap;
0446
0447         break;
0448       }
0449       default:
0450     break;
0451       }
0452     }
0453
0454     // encode layer map
0455     layers = (mvtx_l[0] ? 1U : 0)     // L0 → bit 0
0456              | (mvtx_l[1] ? 2U : 0)   // L1 → bit 1
0457              | (mvtx_l[2] ? 4U : 0);  // L2 → bit 2
0458
0459     // encode state map
0460     states = (mvtx_l_state[0] ? 1U : 0)     // L0 → bit 0
0461              | (mvtx_l_state[1] ? 2U : 0)   // L1 → bit 1
0462              | (mvtx_l_state[2] ? 4U : 0);  // L2 → bit 2
0463     tree->Fill();
0464   }
0465
0466   // clear event vectors for tracks
0467   cluster_vector.clear();
0468   duplicate_cluster_vector.clear();
0469   cluster_vector_INTT.clear();
0470   duplicate_cluster_vector_INTT.clear();
0471
0472   return Fun4AllReturnCodes::EVENT_OK;
0473 }
0474
0475 //____________________________________________________________________________..
0476 int MvtxMatchingEfficiencyWithShapes::EndRun(const int /*runnumber*/)
0477 {
0478   // TFile *root_out = new TFile("MVTX_ME.root","RECREATE");
0479
0480   std::cout << "MvtxMatchingEfficiencyWithShapes::End - Output to " << m_outputFileName << std::endl;
0481
0482   if (PHTFileServer::cd(m_outputFileName))
0483   {
0484     h_status->Write();
0485     h_INTT_time_delta->Write();
0486     tree->Write();
0487   }
0488
0489   std::cout << "MvtxMatchingEfficiencyWithShapes::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0490
0491   return Fun4AllReturnCodes::EVENT_OK;
0492 }
0493
0494 float MvtxMatchingEfficiencyWithShapes::calc_dedx(TrackSeed* tpcseed)
0495 {
0496   std::vector<TrkrDefs::cluskey> clusterKeys;
0497   clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(), tpcseed->end_cluster_keys());
0498
0499   std::vector<float> dedxlist;
0500   for (unsigned long cluster_key : clusterKeys)
0501   {
0502     unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
0503     if (TrkrDefs::getTrkrId(cluster_key) != TrkrDefs::TrkrId::tpcId)
0504     {
0505       continue;
0506     }
0507     TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
0508
0509     float adc = cluster->getAdc();
0510     PHG4TpcCylinderGeom* GeoLayer_local = _geom_container->GetLayerCellGeom(layer_local);
0511     float thick = GeoLayer_local->get_thickness();
0512
0513     float r = GeoLayer_local->get_radius();
0514     float alpha = (r * r) / (2 * r * std::abs(1.0 / tpcseed->get_qOverR()));
0515     float beta = std::atan(tpcseed->get_slope());
0516     float alphacorr = std::cos(alpha);
0517     if (alphacorr < 0 || alphacorr > 4)
0518     {
0519       alphacorr = 4;
0520     }
0521     float betacorr = std::cos(beta);
0522     if (betacorr < 0 || betacorr > 4)
0523     {
0524       betacorr = 4;
0525     }
0526     adc /= thick;
0527     adc *= alphacorr;
0528     adc *= betacorr;
0529     dedxlist.push_back(adc);
0530     sort(dedxlist.begin(), dedxlist.end());
0531   }
0532   int trunc_min = 0;
0533   int trunc_max = (int) dedxlist.size() * 0.7;
0534   float sumdedx = 0;
0535   int ndedx = 0;
0536   for (int j = trunc_min; j <= trunc_max; j++)
0537   {
0538     sumdedx += dedxlist.at(j);
0539     ndedx++;
0540   }
0541   sumdedx /= ndedx;
0542   return sumdedx;
0543 }
0544
0545 std::set<TrkrDefs::cluskey> MvtxMatchingEfficiencyWithShapes::findDuplicates(std::vector<TrkrDefs::cluskey> vec)
0546 {
0547   std::set<TrkrDefs::cluskey> duplicates;
0548   std::sort(vec.begin(), vec.end());
0549   std::set<TrkrDefs::cluskey> distinct(vec.begin(), vec.end());
0550   std::set_difference(vec.begin(), vec.end(), distinct.begin(), distinct.end(),
0551                       std::inserter(duplicates, duplicates.end()));
0552   return duplicates;
0553 }