packages/intt/InttZVertexFinder.cc

0001 #include "InttZVertexFinder.h"
0002 #include "INTTZvtx.h"
0003 #include "InttVertexMapv1.h"
0004 #include "InttVertexv1.h"
0005
0006 /// Fun4All includes
0007 #include <trackbase/ActsGeometry.h>
0008 #include <trackbase/InttDefs.h>
0009 #include <trackbase/TrkrCluster.h>
0010 #include <trackbase/TrkrClusterContainer.h>
0011 #include <trackbase/TrkrHit.h>
0012 #include <trackbase/TrkrHitSet.h>
0013 #include <trackbase/TrkrHitSetContainer.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/getClass.h>
0019
0020 #include <cmath>
0021 #include <filesystem>
0022
0023 /////////////////////////////////////////
0024 // init value for temporal use
0025 std::string sDataType = "data";  // or "MC"
0026 std::string sOutFolderDirectory = "./";
0027 std::pair<double, double> beam_origin = {-0.23436750, 2.5985125};  // note : for run20869
0028 double phi_diff_cut = 1;                                           // 0.265 + 0.269;            // note : for run20869
0029 // std::pair<double, double> DCA_cut             = {0.277 - 0.730, 0.277 + 0.730}; // note : for run20869
0030 std::pair<double, double> DCA_cut = {-3, 3};  // note : for run20869
0031 int N_clu_cut = 10000;
0032 int N_clu_cutl = 4;                 // 20
0033 unsigned int zvtx_cal_require = 3;  // 5; //15;
0034 std::pair<double, double> zvtx_QA_width = {40, 70};
0035 double zvtx_QA_ratio = 0.00006;
0036 bool draw_event_display = false;
0037 bool enable_qa = false;
0038
0039 int clu_sum_adc_cut = 31;
0040 int clu_size_cut = 4;
0041
0042 int data_type = 2;  // note : private gen data cluster
0043
0044 /*
0045 bool GetPhiCheckTag(
0046         std::vector<INTTZvtx::clu_info>& temp_sPH_inner_nocolumn_vec,
0047         std::vector<INTTZvtx::clu_info>& temp_sPH_outer_nocolumn_vec
0048       )
0049 {
0050     int inner_1_check = 0;
0051     int inner_2_check = 0;
0052     int inner_3_check = 0;
0053     int inner_4_check = 0;
0054     for (unsigned int inner_i = 0; inner_i < temp_sPH_inner_nocolumn_vec.size(); inner_i++) {
0055         if (temp_sPH_inner_nocolumn_vec[inner_i].phi >=   0 && temp_sPH_inner_nocolumn_vec[inner_i].phi <  90) inner_1_check = 1;
0056         if (temp_sPH_inner_nocolumn_vec[inner_i].phi >=  90 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 180) inner_2_check = 1;
0057         if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 180 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 270) inner_3_check = 1;
0058         if (temp_sPH_inner_nocolumn_vec[inner_i].phi >= 270 && temp_sPH_inner_nocolumn_vec[inner_i].phi < 360) inner_4_check = 1;
0059
0060         if ( (inner_1_check + inner_2_check + inner_3_check + inner_4_check) == 4 ) break;
0061     }
0062
0063     int outer_1_check = 0;
0064     int outer_2_check = 0;
0065     int outer_3_check = 0;
0066     int outer_4_check = 0;
0067     for (unsigned int outer_i = 0; outer_i < temp_sPH_outer_nocolumn_vec.size(); outer_i++) {
0068         if (temp_sPH_outer_nocolumn_vec[outer_i].phi >=   0 && temp_sPH_outer_nocolumn_vec[outer_i].phi <  90) outer_1_check = 1;
0069         if (temp_sPH_outer_nocolumn_vec[outer_i].phi >=  90 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 180) outer_2_check = 1;
0070         if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 180 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 270) outer_3_check = 1;
0071         if (temp_sPH_outer_nocolumn_vec[outer_i].phi >= 270 && temp_sPH_outer_nocolumn_vec[outer_i].phi < 360) outer_4_check = 1;
0072
0073         if ( (outer_1_check + outer_2_check + outer_3_check + outer_4_check) == 4 ) break;
0074     }
0075
0076     if ( (inner_1_check + inner_2_check + inner_3_check + inner_4_check +
0077           outer_1_check + outer_2_check + outer_3_check + outer_4_check) != 8 ) {return false;}
0078     else { return true; }
0079 }
0080 */
0081
0082 //____________________________________________________________________________..
0083 InttZVertexFinder::InttZVertexFinder(const std::string& name)
0084   : SubsysReco(name)
0085   , m_inttzvtx(new INTTZvtx(sDataType,
0086                             sOutFolderDirectory,
0087                             beam_origin,
0088                             phi_diff_cut,
0089                             DCA_cut,
0090                             N_clu_cutl,
0091                             N_clu_cut,
0092                             zvtx_cal_require,
0093                             zvtx_QA_width,
0094                             draw_event_display,
0095                             enable_qa))
0096 {
0097   std::cout << "InttZVertexFinder::InttZVertexFinder(const std::string &name) Calling ctor" << std::endl;
0098 }
0099
0100 //____________________________________________________________________________..
0101 InttZVertexFinder::~InttZVertexFinder()
0102 {
0103   std::cout << "InttZVertexFinder::~InttZVertexFinder() Calling dtor" << std::endl;
0104   delete m_inttzvtx;
0105 }
0106
0107 //____________________________________________________________________________..
0108 int InttZVertexFinder::Init(PHCompositeNode* /*topNode*/)
0109 {
0110   std::cout << "InttZVertexFinder::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0111   // m_inttzvtx->Init();
0112   //
0113   return Fun4AllReturnCodes::EVENT_OK;
0114 }
0115
0116 //____________________________________________________________________________..
0117 int InttZVertexFinder::InitRun(PHCompositeNode* topNode)
0118 {
0119   std::cout << "InttZVertexFinder::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0120
0121   if (createNodes(topNode) == Fun4AllReturnCodes::ABORTEVENT)
0122   {
0123     return Fun4AllReturnCodes::ABORTEVENT;
0124   }
0125
0126   m_inttzvtx->SetPrintMessageOpt(Verbosity() > 0);
0127
0128   m_inttzvtx->Init();
0129
0130   return Fun4AllReturnCodes::EVENT_OK;
0131 }
0132
0133 int InttZVertexFinder::createNodes(PHCompositeNode* topNode)
0134 {
0135   PHNodeIterator iter(topNode);
0136   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0137   if (!dstNode)
0138   {
0139     std::cout << PHWHERE << "DST Node missing doing nothing" << std::endl;
0140     return Fun4AllReturnCodes::ABORTEVENT;
0141   }
0142
0143   // ---
0144   PHCompositeNode* inttNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "INTT"));
0145   if (!inttNode)
0146   {
0147     inttNode = new PHCompositeNode("INTT");
0148     dstNode->addNode(inttNode);
0149   }
0150
0151   m_inttvertexmap = findNode::getClass<InttVertexMap>(inttNode, "InttVertexMap");
0152   if (!m_inttvertexmap)
0153   {
0154     m_inttvertexmap = new InttVertexMapv1();
0155     PHIODataNode<PHObject>* VertexMapNode = new PHIODataNode<PHObject>(m_inttvertexmap, "InttVertexMap", "PHObject");
0156     inttNode->addNode(VertexMapNode);
0157   }
0158
0159   return Fun4AllReturnCodes::EVENT_OK;
0160 }
0161
0162 //____________________________________________________________________________..
0163 int InttZVertexFinder::process_event(PHCompositeNode* topNode)
0164 {
0165   if (Verbosity() > 1)
0166   {
0167     std::cout << "InttZVertexFinder::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0168   }
0169
0170   //  INTTClu -> EvtInit(event_i);
0171   //  INTTClu -> EvtSetCluGroup();
0172
0173   static int event_i = 0;
0174
0175   ActsGeometry* tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0176   if (!tGeometry)
0177   {
0178     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing." << std::endl;
0179     return Fun4AllReturnCodes::ABORTEVENT;
0180   }
0181
0182   TrkrClusterContainer* clusterMap =
0183       findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0184   if (!clusterMap)
0185   {
0186     std::cout << PHWHERE << "TrkrClusterContainer node is missing." << std::endl;
0187     return Fun4AllReturnCodes::ABORTEVENT;
0188   }
0189
0190   /////////////////////////
0191   // cluster list
0192   //  fill cluster
0193   std::vector<INTTZvtx::clu_info> temp_sPH_inner_nocolumn_vec;
0194   std::vector<INTTZvtx::clu_info> temp_sPH_outer_nocolumn_vec;
0195
0196   std::vector<std::vector<double> > temp_sPH_nocolumn_vec(2);
0197   std::vector<std::vector<double> > temp_sPH_nocolumn_rz_vec(2);
0198
0199   for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++)
0200   {
0201     std::vector<INTTZvtx::clu_info>* p_temp_sPH_nocolumn_vec =
0202         (inttlayer < 2) ? (&temp_sPH_inner_nocolumn_vec)
0203                         : (&temp_sPH_outer_nocolumn_vec);
0204
0205     for (const auto& hitsetkey : clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0206     {
0207       auto range = clusterMap->getClusters(hitsetkey);
0208
0209       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0210       {
0211         const auto cluskey = clusIter->first;
0212         const auto cluster = clusIter->second;
0213
0214         const auto globalPos = tGeometry->getGlobalPosition(cluskey, cluster);
0215         double clu_x = globalPos.x() * 10.;  // convert from "cm" to "mm" unit
0216         double clu_y = globalPos.y() * 10.;  // convert from "cm" to "mm" unit
0217         double clu_z = globalPos.z() * 10.;  // convert from "cm" to "mm" unit
0218
0219         double clu_phi = (clu_y < 0) ? atan2(clu_y, clu_x) * (180. / M_PI) + 360
0220                                      : atan2(clu_y, clu_x) * (180. / M_PI);
0221
0222         double clu_radius = sqrt(pow(clu_x, 2) + pow(clu_y, 2));
0223
0224         // int size = cluster->getSize();
0225         int adc = cluster->getAdc();
0226
0227         if (adc < 40)
0228         {
0229           continue;
0230         }
0231
0232         p_temp_sPH_nocolumn_vec->push_back({-1,
0233                                             -1,
0234                                             (int) cluster->getAdc(),
0235                                             (int) cluster->getAdc(),
0236                                             (int) cluster->getPhiSize(),
0237                                             clu_x,
0238                                             clu_y,
0239                                             clu_z,
0240                                             (int) inttlayer + 3,  // note : should be 3 or 4, for the inner layer, this is for the mega cluster search
0241                                             clu_phi});
0242
0243         temp_sPH_nocolumn_vec[0].push_back(clu_x);
0244         temp_sPH_nocolumn_vec[1].push_back(clu_y);
0245
0246         temp_sPH_nocolumn_rz_vec[0].push_back(clu_z);
0247         temp_sPH_nocolumn_rz_vec[1].push_back((clu_phi > 180) ? clu_radius * -1 : clu_radius);
0248
0249         // std::cout<<"cluster pos: ";
0250         // std::cout<<globalPos.x()<<" ";
0251         // std::cout<<globalPos.y()<<" ";
0252         // std::cout<<globalPos.z()<<std::endl;
0253       }
0254     }
0255   }
0256
0257   ////////////////////////
0258
0259   int NvtxMC = 1;
0260   double TrigZvtxMC = 0.;
0261   Long64_t bcofull = 0;
0262
0263   bool status = m_inttzvtx->ProcessEvt(
0264       event_i,
0265       temp_sPH_inner_nocolumn_vec,
0266       temp_sPH_outer_nocolumn_vec,
0267       temp_sPH_nocolumn_vec,
0268       temp_sPH_nocolumn_rz_vec,
0269       NvtxMC,
0270       TrigZvtxMC,
0271       true,  // GetPhiCheckTag(temp_sPH_inner_nocolumn_vec, temp_sPH_outer_nocolumn_vec),
0272       bcofull,
0273       5  // centrality bin, note : no bco_full for MC
0274   );
0275
0276   std::cout << "InttZVertex:process_evt status = " << (status ? "good" : "failed") << std::endl;
0277
0278   std::vector<double> vtxout = m_inttzvtx->GetEvtZPeak();              // mm unit
0279   std::pair<double, double> beamorigin = m_inttzvtx->GetBeamOrigin();  // mm unit
0280
0281   std::vector<double> final_vtx{beamorigin.first * 0.1,
0282                                 beamorigin.second * 0.1,
0283                                 vtxout[1] * 0.1};  // mm -> cm by 0.1
0284
0285   INTTZvtx::ZvtxInfo& zvtxinfo = m_inttzvtx->GetZvtxInfo();
0286
0287   //--if(m_inttvertexmap->size()>0){
0288   //--  std::cout<<"added to existing vertex"<<std::endl;
0289   //--  for(auto vtx = m_inttvertexmap->begin(); vtx!=m_inttvertexmap->end(); ++vtx){
0290   //--    vtx->second->set_x(final_vtx[0]);
0291   //--    vtx->second->set_y(final_vtx[1]);
0292   //--    vtx->second->set_z(final_vtx[2]);
0293   //--  }
0294   //--}
0295   //--else {
0296   //--std::cout<<"new vertex created"<<std::endl;
0297   auto vertex = std::make_unique<InttVertexv1>();
0298   vertex->set_x(final_vtx[0]);
0299   vertex->set_y(final_vtx[1]);
0300   vertex->set_z(final_vtx[2]);
0301
0302   vertex->set_chi2ndf(zvtxinfo.chi2ndf);
0303   vertex->set_width(zvtxinfo.width);
0304   vertex->set_good(zvtxinfo.good);
0305   vertex->set_nclus(zvtxinfo.nclus);
0306   vertex->set_ntracklet(zvtxinfo.ntracklets);
0307   vertex->set_ngroup(zvtxinfo.ngroup);
0308   vertex->set_peakratio(zvtxinfo.peakratio);
0309   vertex->set_peakwidth(zvtxinfo.peakwidth);
0310
0311   // vertex->set_z_err(); // no value asigned yet
0312   m_inttvertexmap->insert(vertex.release());
0313   //}
0314
0315   if (Verbosity() > 0)
0316   {
0317     std::cout << "vecsize in : " << temp_sPH_inner_nocolumn_vec.size()
0318               << ", out: " << temp_sPH_outer_nocolumn_vec.size()
0319               << ", zvtx : " << vtxout[1] << " +- " << vtxout[2]
0320               << " sts : " << vtxout[0]
0321               << std::endl;
0322   }
0323
0324   m_inttzvtx->ClearEvt();
0325
0326   event_i++;
0327
0328   return Fun4AllReturnCodes::EVENT_OK;
0329 }
0330
0331 //____________________________________________________________________________..
0332 int InttZVertexFinder::End(PHCompositeNode* /*topNode*/)
0333 {
0334   if (Verbosity())
0335   {
0336     std::cout << "InttZVertexFinder::End(PHCompositeNode *topNode) " << std::endl;
0337   }
0338
0339   m_inttzvtx->PrintPlots();
0340   m_inttzvtx->EndRun();
0341
0342   return Fun4AllReturnCodes::EVENT_OK;
0343 }
0344
0345 //____________________________________________________________________________..
0346 void InttZVertexFinder::Print(const std::string& what) const
0347 {
0348   std::cout << "InttZVertexFinder::Print(const std::string &what) const Printing info for " << what << std::endl;
0349 }
0350
0351 void InttZVertexFinder::SetBeamCenter(const double beamx, const double beamy)
0352 {
0353   if (m_inttzvtx != nullptr)
0354   {
0355     m_inttzvtx->SetBeamOrigin(beamx * 10., beamy * 10.);  // convert to cm to mm unit
0356   }
0357 }
0358
0359 void InttZVertexFinder::SetOutDirectory(const std::string& outDirectory)
0360 {
0361   if (m_inttzvtx != nullptr)
0362   {
0363     m_inttzvtx->SetOutDirectory(outDirectory);  // convert to cm to mm unit
0364   }
0365 }
0366
0367 void InttZVertexFinder::EnableQA(const bool enableQA)
0368 {
0369   if (m_inttzvtx != nullptr)
0370   {
0371     m_inttzvtx->EnableQA(enableQA);
0372   }
0373 }
0374
0375 void InttZVertexFinder::EnableEventDisplay(const bool enableEvtDisp)
0376 {
0377   if (m_inttzvtx != nullptr)
0378   {
0379     m_inttzvtx->EnableEventDisplay(enableEvtDisp);
0380   }
0381 }