packages/intt/InttXYVertexFinder.cc

0001 #include "InttXYVertexFinder.h"
0002 #include "INTTXYvtx.h"
0003 #include "InttVertexMapv1.h"
0004 #include "InttVertexv1.h"
0005
0006 #include <trackbase/ActsGeometry.h>
0007 #include <trackbase/InttDefs.h>
0008 #include <trackbase/TrkrCluster.h>
0009 #include <trackbase/TrkrClusterContainer.h>
0010 #include <trackbase/TrkrHit.h>
0011 #include <trackbase/TrkrHitSet.h>
0012 #include <trackbase/TrkrHitSetContainer.h>
0013
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/getClass.h>
0018
0019 #include <cmath>
0020 #include <string>
0021
0022 std::string xyinit_sDataType = "data";  // or "MC"
0023 std::string xyinit_input_directory = "./";
0024 std::string xyinit_file_name = "MC_ZF_xyvtx";
0025 std::string xyinit_out_folder_directory = xyinit_input_directory + "/SemiFinal_RunXY_" + xyinit_file_name;
0026
0027 // std::pair<double, double> xyinit_beam_origin = {-0.015, 0.012};
0028
0029 std::pair<double, double> xyinit_beam_origin = {0.0, 0.0};
0030 std::pair<double, double> xyinit_DCA_cut = {-1, 1};  // note : used
0031 double xyinit_phi_diff_cut = 0.35;                   // note : used
0032 int xyinit_clu_sum_adc_cut = -1;                     // note : used
0033 int xyinit_clu_size_cut = 4;                         // note : used
0034 int xyinit_N_clu_cut = 350;                          // note : used
0035 int xyinit_N_clu_cutl = 20;                          // note : used
0036
0037 int xyinit_zvtx_cal_require = 15;
0038 double xyinit_peek = 3.324;
0039
0040 bool GetPhiCheckTag(
0041     std::vector<INTTXYvtx::clu_info>& temp_sPH_inner_nocolumn_vec,
0042     std::vector<INTTXYvtx::clu_info>& temp_sPH_outer_nocolumn_vec)
0043 {
0044   int inner_1_check = 0;
0045   int inner_2_check = 0;
0046   int inner_3_check = 0;
0047   int inner_4_check = 0;
0048   // for (unsigned int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++)
0049   for (const auto& clu : temp_sPH_inner_nocolumn_vec)
0050   {
0051     if (clu.phi >= 0 && clu.phi < 90)
0052     {
0053       inner_1_check = 1;
0054     }
0055     if (clu.phi >= 90 && clu.phi < 180)
0056     {
0057       inner_2_check = 1;
0058     }
0059     if (clu.phi >= 180 && clu.phi < 270)
0060     {
0061       inner_3_check = 1;
0062     }
0063     if (clu.phi >= 270 && clu.phi < 360)
0064     {
0065       inner_4_check = 1;
0066     }
0067
0068     if ((inner_1_check + inner_2_check + inner_3_check + inner_4_check) == 4)
0069     {
0070       break;
0071     }
0072   }
0073
0074   int outer_1_check = 0;
0075   int outer_2_check = 0;
0076   int outer_3_check = 0;
0077   int outer_4_check = 0;
0078   // for (unsigned int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++)
0079   for (const auto& clu : temp_sPH_outer_nocolumn_vec)
0080   {
0081     if (clu.phi >= 0 && clu.phi < 90)
0082     {
0083       outer_1_check = 1;
0084     }
0085     if (clu.phi >= 90 && clu.phi < 180)
0086     {
0087       outer_2_check = 1;
0088     }
0089     if (clu.phi >= 180 && clu.phi < 270)
0090     {
0091       outer_3_check = 1;
0092     }
0093     if (clu.phi >= 270 && clu.phi < 360)
0094     {
0095       outer_4_check = 1;
0096     }
0097
0098     if ((outer_1_check + outer_2_check + outer_3_check + outer_4_check) == 4)
0099     {
0100       break;
0101     }
0102   }
0103
0104   if ((inner_1_check + inner_2_check + inner_3_check + inner_4_check +
0105        outer_1_check + outer_2_check + outer_3_check + outer_4_check) != 8)
0106   {
0107     return false;
0108   }
0109   else
0110   {
0111     return true;
0112   }
0113 }
0114
0115 //____________________________________________________________________________..
0116 InttXYVertexFinder::InttXYVertexFinder(const std::string& name)
0117   : SubsysReco(name)
0118   , m_inttxyvtx(new INTTXYvtx(
0119         xyinit_sDataType, xyinit_out_folder_directory, xyinit_beam_origin, xyinit_phi_diff_cut, xyinit_DCA_cut, xyinit_N_clu_cutl, xyinit_N_clu_cut, xyinit_peek))
0120 {
0121   std::cout << "InttXYVertexFinder::InttXYVertexFinder(const std::string &name) Calling ctor" << std::endl;
0122 }
0123
0124 //____________________________________________________________________________..
0125 InttXYVertexFinder::~InttXYVertexFinder()
0126 {
0127   std::cout << "InttXYVertexFinder::~InttXYVertexFinder() Calling dtor" << std::endl;
0128   delete m_inttxyvtx;
0129 }
0130
0131 //____________________________________________________________________________..
0132 int InttXYVertexFinder::Init(PHCompositeNode* /*topNode*/)
0133 {
0134   std::cout << "InttXYVertexFinder::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0135   return Fun4AllReturnCodes::EVENT_OK;
0136 }
0137
0138 //____________________________________________________________________________..
0139 int InttXYVertexFinder::InitRun(PHCompositeNode* topNode)
0140 {
0141   std::cout << "InttXYVertexFinder::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0142
0143   // should be reset before init
0144   m_inttxyvtx->PrintMessageOpt(Verbosity() > 5);
0145   m_inttxyvtx->Init();
0146
0147   if (createNodes(topNode) == Fun4AllReturnCodes::ABORTEVENT)
0148   {
0149     return Fun4AllReturnCodes::ABORTEVENT;
0150   }
0151
0152   return Fun4AllReturnCodes::EVENT_OK;
0153 }
0154
0155 //____________________________________________________________________________..
0156 int InttXYVertexFinder::process_event(PHCompositeNode* topNode)
0157 {
0158   if (Verbosity() > 1)
0159   {
0160     std::cout << "InttXYVertexFinder::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0161   }
0162
0163   static int event_i = 0;
0164
0165   ActsGeometry* tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0166   if (!tGeometry)
0167   {
0168     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing." << std::endl;
0169     return Fun4AllReturnCodes::ABORTEVENT;
0170   }
0171
0172   TrkrClusterContainer* clusterMap =
0173       findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0174   if (!clusterMap)
0175   {
0176     std::cout << PHWHERE << "TrkrClusterContainer node is missing." << std::endl;
0177     return Fun4AllReturnCodes::ABORTEVENT;
0178   }
0179
0180   /////////////////////////
0181   // cluster list
0182   //  fill cluster
0183   std::vector<INTTXYvtx::clu_info> temp_sPH_inner_nocolumn_vec;
0184   std::vector<INTTXYvtx::clu_info> temp_sPH_outer_nocolumn_vec;
0185
0186   std::vector<std::vector<double>> temp_sPH_nocolumn_vec(2);
0187   std::vector<std::vector<double>> temp_sPH_nocolumn_rz_vec(2);
0188
0189   for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++)
0190   {
0191     std::vector<INTTXYvtx::clu_info>* p_temp_sPH_nocolumn_vec =
0192         (inttlayer < 2) ? (&temp_sPH_inner_nocolumn_vec)
0193                         : (&temp_sPH_outer_nocolumn_vec);
0194
0195     for (const auto& hitsetkey : clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0196     {
0197       auto range = clusterMap->getClusters(hitsetkey);
0198
0199       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0200       {
0201         const auto cluskey = clusIter->first;
0202         const auto cluster = clusIter->second;
0203
0204         const auto globalPos = tGeometry->getGlobalPosition(cluskey, cluster);
0205         double clu_x = globalPos.x() * 10.;  // convert to "mm" unit
0206         double clu_y = globalPos.y() * 10.;  // convert to "mm" unit
0207         double clu_z = globalPos.z() * 10.;  // convert to "mm" unit
0208
0209         double clu_phi = (clu_y < 0) ? atan2(clu_y, clu_x) * (180. / M_PI) + 360
0210                                      : atan2(clu_y, clu_x) * (180. / M_PI);
0211
0212         double clu_radius = sqrt(pow(clu_x, 2) + pow(clu_y, 2));
0213
0214         p_temp_sPH_nocolumn_vec->push_back({-1,
0215                                             -1,
0216                                             (int) cluster->getAdc(),
0217                                             (int) cluster->getAdc(),
0218                                             (int) cluster->getPhiSize(),
0219                                             clu_x,                // change to "mm"
0220                                             clu_y,                // change to "mm"
0221                                             clu_z,                // change to "mm"
0222                                             (int) inttlayer + 3,  // note : should be 3 or 4, for the inner layer, this is for the mega cluster search
0223                                             clu_phi});
0224
0225         temp_sPH_nocolumn_vec[0].push_back(clu_x);
0226         temp_sPH_nocolumn_vec[1].push_back(clu_y);
0227
0228         temp_sPH_nocolumn_rz_vec[0].push_back(clu_z);
0229         temp_sPH_nocolumn_rz_vec[1].push_back((clu_phi > 180) ? clu_radius * -1 : clu_radius);
0230       }
0231     }
0232   }
0233
0234   ////////////////////////
0235   //
0236   int NvtxMC = 1;
0237   double TrigZvtxMC = 0.;
0238   Long64_t bcofull = -1;
0239
0240   m_inttxyvtx->ProcessEvt(
0241       event_i,
0242       temp_sPH_inner_nocolumn_vec,
0243       temp_sPH_outer_nocolumn_vec,
0244       temp_sPH_nocolumn_vec,
0245       temp_sPH_nocolumn_rz_vec,
0246       NvtxMC,
0247       TrigZvtxMC,
0248       GetPhiCheckTag(temp_sPH_inner_nocolumn_vec, temp_sPH_outer_nocolumn_vec),
0249       bcofull  // note : no bco_full for MC
0250   );
0251
0252   m_inttxyvtx->ClearEvt();
0253
0254   ///////////////////////////////////////////////////
0255   // calculate XY vertex
0256   if ((event_i % m_period) == 0)
0257   {
0258     // quadorant method
0259     std::vector<std::pair<double, double>> out_vtx = m_inttxyvtx->MacroVTXSquare(4, 10);
0260     m_vertex_quad[0] = out_vtx[0].first;
0261     m_vertex_quad[1] = out_vtx[0].second;
0262
0263     if (Verbosity() > 1)
0264     {
0265       std::cout << " " << std::endl;
0266       std::cout << "The best vertex throughout the scan: " << out_vtx[0].first << " " << out_vtx[0].second << std::endl;
0267       std::cout << "The origin during that scan: " << out_vtx[1].first << " " << out_vtx[1].second << std::endl;
0268       std::cout << "Fit error, DCA and angle diff: " << out_vtx[2].first << " " << out_vtx[2].second << std::endl;
0269       std::cout << "fit pol0 pos Y, DCA and angle diff: " << out_vtx[3].first << " " << out_vtx[3].second << std::endl;
0270     }
0271
0272     // line filled method
0273     std::vector<std::pair<double, double>> out_vtx_line = m_inttxyvtx->FillLine_FindVertex(
0274         {(out_vtx[0].first + out_vtx[1].first) / 2.,
0275          (out_vtx[0].second + out_vtx[1].second) / 2.},
0276         0.001);
0277
0278     m_vertex_line[0] = out_vtx_line[0].first;
0279     m_vertex_line[1] = out_vtx_line[0].second;
0280
0281     if (Verbosity() > 1)
0282     {
0283       std::cout << " " << std::endl;
0284       std::cout << "By fill-line method," << std::endl;
0285       std::cout << "Reco Run Vertex XY: " << out_vtx_line[0].first << " " << out_vtx_line[0].second << std::endl;
0286       std::cout << "Reco Run Vertex XY Error: " << out_vtx_line[1].first << " " << out_vtx_line[1].second << std::endl;
0287     }
0288
0289     ///////////////////////////////////////////////////
0290     // convert xy-vertex from "mm" to "cm"  unit
0291     m_vertex_quad[0] *= 0.1;
0292     m_vertex_quad[1] *= 0.1;
0293     m_vertex_line[0] *= 0.1;
0294     m_vertex_line[1] *= 0.1;
0295   }
0296
0297   auto vertex_quad = std::make_unique<InttVertexv1>();
0298   vertex_quad->set_x(m_vertex_quad[0]);
0299   vertex_quad->set_y(m_vertex_quad[1]);
0300   vertex_quad->set_z(m_vertex_quad[2]);
0301   // vertex_quad->set_id(0); // automatically set in insert
0302   //  vertex->set_z_err(); // no value asigned yet
0303   m_inttvertexmap->insert(vertex_quad.release());
0304
0305   auto vertex_line = std::make_unique<InttVertexv1>();
0306   vertex_line->set_x(m_vertex_line[0]);
0307   vertex_line->set_y(m_vertex_line[1]);
0308   vertex_line->set_z(m_vertex_line[2]);
0309   // vertex_line->set_id(1); // automatically set in insert
0310   //  vertex->set_z_err(); // no value asigned yet
0311   m_inttvertexmap->insert(vertex_line.release());
0312
0313   if (
0314       Verbosity() > 0 &&
0315       (event_i % m_period) == 0)
0316   {
0317     std::cout << event_i << ", vtx_xy_quad: " << m_vertex_quad[0] << " " << m_vertex_quad[1] << " " << m_vertex_quad[2]
0318               << ", vtx_xy_line: " << m_vertex_line[0] << " " << m_vertex_line[1] << " " << m_vertex_line[2]
0319               << std::endl;
0320   }
0321   event_i++;
0322
0323   return Fun4AllReturnCodes::EVENT_OK;
0324 }
0325
0326 //____________________________________________________________________________..
0327 int InttXYVertexFinder::End(PHCompositeNode* /*topNode*/)
0328 {
0329   m_inttxyvtx->PrintPlots();
0330   m_inttxyvtx->EndRun();
0331
0332   return Fun4AllReturnCodes::EVENT_OK;
0333 }
0334
0335 //____________________________________________________________________________..
0336 void InttXYVertexFinder::Print(const std::string& what) const
0337 {
0338   std::cout << "InttXYVertexFinder::Print(const std::string &what) const Printing info for " << what << std::endl;
0339 }
0340
0341 int InttXYVertexFinder::createNodes(PHCompositeNode* topNode)
0342 {
0343   PHNodeIterator iter(topNode);
0344   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0345   if (!dstNode)
0346   {
0347     std::cout << PHWHERE << "DST Node missing doing nothing" << std::endl;
0348     return Fun4AllReturnCodes::ABORTEVENT;
0349   }
0350
0351   // ---
0352   PHCompositeNode* inttNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "INTT"));
0353   if (!inttNode)
0354   {
0355     inttNode = new PHCompositeNode("INTT");
0356     dstNode->addNode(inttNode);
0357   }
0358
0359   m_inttvertexmap = findNode::getClass<InttVertexMap>(inttNode, "InttVertexMap");
0360   if (!m_inttvertexmap)
0361   {
0362     m_inttvertexmap = new InttVertexMapv1();
0363     PHIODataNode<PHObject>* VertexMapNode = new PHIODataNode<PHObject>(m_inttvertexmap, "InttVertexMap", "PHObject");
0364     inttNode->addNode(VertexMapNode);
0365   }
0366
0367   return Fun4AllReturnCodes::EVENT_OK;
0368 }
0369
0370 void InttXYVertexFinder::SetBeamCenter(const double beamx, const double beamy)
0371 {
0372   if (m_inttxyvtx != nullptr)
0373   {
0374     m_inttxyvtx->SetBeamOrigin(beamx * 10., beamy * 10.);  // converted from "cm" to "mm"
0375   }
0376 }
0377
0378 void InttXYVertexFinder::SetSaveHisto(const bool savehist)
0379 {
0380   if (m_inttxyvtx != nullptr)
0381   {
0382     m_inttxyvtx->SetSaveHisto(savehist);
0383   }
0384 }
0385
0386 void InttXYVertexFinder::EnableDrawHisto(const bool enable)
0387 {
0388   if (m_inttxyvtx != nullptr)
0389   {
0390     m_inttxyvtx->EnableDrawHisto(enable);
0391   }
0392 }
0393
0394 void InttXYVertexFinder::EnableQA(const bool enable)
0395 {
0396   if (m_inttxyvtx != nullptr)
0397   {
0398     m_inttxyvtx->EnableQA(enable);
0399   }
0400 }