SiCalo/SiliconSeedAna/ChargeRecalc.h

0001 #ifndef __CHARGE_RECALC_H__
0002 #define __CHARGE_RECALC_H__
0003
0004 #include <fun4all/SubsysReco.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/getClass.h>
0009
0010 #include <trackbase/ActsGeometry.h>
0011 #include <trackbase/TrkrClusterContainer.h>
0012
0013 #include <trackbase_historic/SvtxTrackMap.h>
0014 #include <trackbase_historic/SvtxTrack.h>
0015 #include <trackbase_historic/TrackSeed.h> // explicit (avoid relying on indirect includes)
0016
0017 #include <globalvertex/SvtxVertexMap.h>
0018 #include <globalvertex/SvtxVertex.h>
0019
0020 #include <cmath>
0021 #include <limits>
0022 #include <iostream>
0023
0024 static inline float WrapPhi(float x)
0025 {
0026     while (x > M_PI)
0027         x -= 2.0f * M_PI;
0028     while (x < -M_PI)
0029         x += 2.0f * M_PI;
0030     return x;
0031 }
0032
0033 static inline float DeltaPhi(float a, float b)
0034 {
0035     return WrapPhi(a - b);
0036 }
0037
0038 class ChargeRecalc : public SubsysReco
0039 {
0040 public:
0041     ChargeRecalc(const std::string &name = "ChargeRecalc")
0042         : SubsysReco(name) {}
0043     void setVerbosity(int v) { m_verb = v; }
0044
0045     int process_event(PHCompositeNode *topNode) override
0046     {
0047         auto clustermap = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
0048         auto geometry = findNode::getClass<ActsGeometry>(topNode, m_actsgeometryName);
0049         auto vertexmap = findNode::getClass<SvtxVertexMap>(topNode, m_vertexMapName);
0050         auto trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0051         if (!trackmap)
0052         {
0053             std::cout << PHWHERE << "Missing trackmap, can't continue" << std::endl;
0054             return Fun4AllReturnCodes::EVENT_OK;
0055         }
0056         if (!clustermap)
0057         {
0058             std::cout << PHWHERE << "Missing clustermap, can't continue" << std::endl;
0059             return Fun4AllReturnCodes::EVENT_OK;
0060         }
0061         if (!geometry)
0062         {
0063             std::cout << PHWHERE << "Missing geometry, can't continue" << std::endl;
0064             return Fun4AllReturnCodes::EVENT_OK;
0065         }
0066
0067         // if (m_verb > 0)
0068         // {
0069         //     float vx = NAN, vy = NAN, vz = NAN;
0070         //     int vc = 999;
0071
0072              if (vertexmap && !vertexmap->empty())
0073              {
0074                std::cout<<"Ntrk: "<<trackmap->size()<<", Nvtx: "<<vertexmap->size()<<"  "<<std::endl;
0075                int ivtx=0;
0076                for (const auto &[k, v] : *vertexmap)
0077                {
0078                   if (v){
0079                    std::cout<<"vtx : "<< ivtx  <<" "
0080                        << v->get_id() <<" "
0081                        << v->get_x() <<" "
0082                        << v->get_y() <<" "
0083                        << v->get_z() <<" "
0084                        << v->get_beam_crossing() <<std::endl;
0085                   }
0086                   ivtx++;
0087                }
0088         //         SvtxVertex *vtx = vertexmap->begin()->second;
0089         //         if (vtx)
0090         //         {
0091         //             vx = vtx->get_x();
0092         //             vy = vtx->get_y();
0093         //             vz = vtx->get_z();
0094         //             vc = vtx->get_beam_crossing();
0095         //         }
0096              }
0097
0098         //     std::cout << "[ChargeRecalc] vtx(first): x=" << vx
0099         //               << " y=" << vy << " z=" << vz
0100         //               << " crossing=" << vc
0101         //               << std::endl;
0102         // }
0103
0104         int n_changed = 0;
0105         int n_skipped = 0;
0106         int n_total = 0;
0107
0108         for (auto &iter : *trackmap)
0109         {
0110             SvtxTrack *track = iter.second;
0111             if (!track)
0112                 continue;
0113             n_total++;
0114
0115             int vertexid = track->get_vertex_id();
0116             std::cout<<"vertexID : "<<vertexid<<std::endl;
0117
0118             TrackSeed *si_seed = track->get_silicon_seed();
0119             if (!si_seed)
0120             {
0121                 n_skipped++;
0122                 continue;
0123             }
0124
0125 //--
0126 //--          float xvtx = 0.0f;
0127 //--          float yvtx = 0.0f;
0128 //--
0129 //--          bool found_vtx = false;
0130 //--          if (vertexmap && !vertexmap->empty())
0131 //--          {
0132 //--              SvtxVertex *vtx0 = nullptr;
0133 //--
0134 //--              // crossing==0 has high priority
0135 //--              for (const auto &[k, v] : *vertexmap)
0136 //--              {
0137 //--                  if (v && v->get_beam_crossing() == 0)
0138 //--                  {
0139 //--                      vtx0 = v;
0140 //--                      break;
0141 //--                  }
0142 //--              }
0143 //--
0144 //--              if (!vtx0)
0145 //--                  vtx0 = vertexmap->begin()->second;
0146 //--
0147 //--              if (vtx0)
0148 //--              {
0149 //--                  xvtx = vtx0->get_x();
0150 //--                  yvtx = vtx0->get_y();
0151 //--                  //--float zvtx = vtx0->get_z();
0152 //--                  //--std::cout<<"Vtx : "<<xvtx<<" "<<yvtx<<" "<<zvtx<<std::endl;
0153 //--                  found_vtx = true;
0154 //--              }
0155 //--          }
0156 //--
0157 //--          if (m_verb > 0 && !found_vtx)
0158 //--          {
0159 //--              std::cout << "[ChargeRecalc] WARNING: no usable vertex (vertexmap="
0160 //--                        << (void *)vertexmap
0161 //--                        << ", size=" << (vertexmap ? (int)vertexmap->size() : -1)
0162 //--                        << ")" << std::endl;
0163 //--              //--n_skipped++;
0164 //--              //--continue;
0165 //--          }
0166 //--
0167
0168             // find innermost/outermost cluster phi from vertex
0169 //--            float rmin = std::numeric_limits<float>::max();
0170 //--            float rmax = -1.0f;
0171             float phi_rmin = 0.0f;
0172             float phi_rmax = 0.0f;
0173
0174             std::vector< Acts::Vector3 > vClus;
0175             for (auto key_iter = si_seed->begin_cluster_keys(); key_iter != si_seed->end_cluster_keys(); ++key_iter)
0176             {
0177                 const auto &cluster_key = *key_iter;
0178                 auto trkrCluster = clustermap->findCluster(cluster_key);
0179                 if (!trkrCluster)
0180                     continue;
0181
0182                 int  layer =  TrkrDefs::getLayer(cluster_key);
0183                 auto global = geometry->getGlobalPosition(cluster_key, trkrCluster);
0184
0185                 vClus.push_back(global);
0186
0187                 float clus_r = sqrt(global[0]*global[0] + global[1]*global[1]);
0188                 std::cout<<"clus : "<<layer<<" "<<global[0]<<" "<<global[1]<<" "<<clus_r<<std::endl;
0189
0190                 // original algorithm
0191                 //--const float dx = global[0] - xvtx;
0192                 //--const float dy = global[1] - yvtx;
0193                 //--const float r = sqrt(dx * dx + dy * dy);
0194                 //--const float phi = std::atan2(dy, dx);
0195
0196                 //--if (r < rmin)
0197                 //--{
0198                 //--    rmin = r;
0199                 //--    phi_rmin = phi;
0200                 //--}
0201                 //--if (r > rmax)
0202                 //--{
0203                 //--    rmax = r;
0204                 //--    phi_rmax = phi;
0205                 //--}
0206             }
0207
0208             // original algorithm
0209             //--if (rmax < 0.0f)
0210             //--{
0211             //--    n_skipped++;
0212             //--    continue;
0213             //--}
0214
0215             if(vClus.size()<3) { std::cout<<"Ncluster is "<<vClus.size()<<" <= 3. skip"<<std::endl;  continue; }
0216
0217             phi_rmin = atan2(vClus[1][1]-vClus[0][1], vClus[1][0]-vClus[0][0]);
0218             phi_rmax = atan2(vClus[vClus.size()-1][1]-vClus[0][1], vClus[vClus.size()-1][0]-vClus[0][0]);
0219
0220             const float dphi = DeltaPhi(phi_rmax, phi_rmin);
0221             const int q_new = (dphi > 0.0f) ? -1 : +1; // if official is flipped, swap -1/+1
0222             const int q_old = track->get_charge();
0223
0224             if (q_new != q_old)
0225             {
0226                 if (m_verb > 2)
0227                 {
0228                     std::cout << "[ChargeRecalc] trkid=" << track->get_id()
0229                               << " before=" << track->get_charge();
0230                 }
0231                 track->set_charge(q_new);
0232                 si_seed->set_qOverR(q_new * fabs(si_seed->get_qOverR()));
0233                 if (m_verb > 2)
0234                 {
0235                     std::cout << " after=" << track->get_charge() << std::endl;
0236                 }
0237                 n_changed++;
0238             }
0239         }
0240         if (m_verb > 1)
0241         {
0242             std::cout << "[ChargeRecalc] total=" << n_total
0243                       << " changed=" << n_changed
0244                       << " skipped=" << n_skipped
0245                       << std::endl;
0246         }
0247         return Fun4AllReturnCodes::EVENT_OK;
0248     }
0249
0250 private:
0251     // node names (match SiliconSeedsAna default)
0252     std::string m_clusterContainerName = "TRKR_CLUSTER";
0253     std::string m_actsgeometryName = "ActsGeometry";
0254     std::string m_trackMapName = "SvtxTrackMap";
0255     std::string m_vertexMapName = "SvtxVertexMap";
0256
0257     int m_verb = 0;
0258 };
0259
0260 #endif