Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:19:58 

0001 #include "PhotonClusterBuilder.h"
0002 
0003 #include <calobase/PhotonClusterv1.h>
0004 #include <calobase/RawCluster.h>
0005 #include <calobase/RawClusterContainer.h>
0006 #include <calobase/RawClusterUtility.h>
0007 #include <calobase/RawTowerGeomContainer.h>
0008 #include <calobase/TowerInfoContainer.h>
0009 #include <calobase/TowerInfoDefs.h>
0010 
0011 // Tower stuff
0012 #include <calobase/RawTowerGeom.h>
0013 #include <calobase/TowerInfo.h>
0014 
0015 // for the vertex
0016 #include <globalvertex/GlobalVertex.h>
0017 #include <globalvertex/GlobalVertexMap.h>
0018 #include <globalvertex/MbdVertex.h>
0019 #include <globalvertex/MbdVertexMap.h>
0020 
0021 #include <phool/PHCompositeNode.h>
0022 #include <phool/PHIODataNode.h>
0023 #include <phool/PHNodeIterator.h>
0024 #include <phool/PHObject.h>
0025 #include <phool/getClass.h>
0026 
0027 #include <fun4all/Fun4AllReturnCodes.h>
0028 
0029 #include <TMVA/RBDT.hxx>
0030 
0031 #include <algorithm>
0032 #include <cmath>
0033 #include <iostream>
0034 #include <set>
0035 #include <stdexcept>
0036 
0037 namespace
0038 {
0039   // Helper function to shift tower indices for wrapping in phi
0040   void shift_tower_index(int& ieta, int& iphi, int etadiv, int phidiv)
0041   {
0042     while (iphi < 0)
0043     {
0044       iphi += phidiv;
0045     }
0046     while (iphi >= phidiv)
0047     {
0048       iphi -= phidiv;
0049     }
0050     if (ieta < 0 || ieta >= etadiv)
0051     {
0052       ieta = -1;  // invalid
0053     }
0054   }
0055 }  // namespace
0056 
0057 PhotonClusterBuilder::PhotonClusterBuilder(const std::string& name)
0058   : SubsysReco(name)
0059 {
0060 }
0061 
0062 int PhotonClusterBuilder::InitRun(PHCompositeNode* topNode)
0063 {
0064   // BDT
0065   if (m_do_bdt)
0066   {
0067     m_bdt = std::make_unique<TMVA::Experimental::RBDT>("myBDT", m_bdt_model_file);
0068   }
0069 
0070   // locate input raw cluster container
0071   m_rawclusters = findNode::getClass<RawClusterContainer>(topNode, m_input_cluster_node);
0072   if (!m_rawclusters)
0073   {
0074     std::cerr << Name() << ": could not find RawClusterContainer node '" << m_input_cluster_node << "'" << std::endl;
0075     return Fun4AllReturnCodes::ABORTRUN;
0076   }
0077 
0078   m_emc_tower_container = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0079   if (!m_emc_tower_container)
0080   {
0081     std::cerr << Name() << ": could not find TowerInfoContainer node 'TOWERINFO_CALIB_CEMC'" << std::endl;
0082     return Fun4AllReturnCodes::ABORTRUN;
0083   }
0084 
0085   m_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0086   if (!m_geomEM)
0087   {
0088     std::cerr << Name() << ": could not find RawTowerGeomContainer node 'TOWERGEOM_CEMC'" << std::endl;
0089     return Fun4AllReturnCodes::ABORTRUN;
0090   }
0091 
0092   m_ihcal_tower_container = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0093   if (!m_ihcal_tower_container)
0094   {
0095     std::cerr << Name() << ": could not find TowerInfoContainer node 'TOWERINFO_CALIB_HCALIN'" << std::endl;
0096     return Fun4AllReturnCodes::ABORTRUN;
0097   }
0098 
0099   m_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0100   if (!m_geomIH)
0101   {
0102     std::cerr << Name() << ": could not find RawTowerGeomContainer node 'TOWERGEOM_HCALIN'" << std::endl;
0103     return Fun4AllReturnCodes::ABORTRUN;
0104   }
0105 
0106   m_ohcal_tower_container = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0107   if (!m_ohcal_tower_container)
0108   {
0109     std::cerr << Name() << ": could not find TowerInfoContainer node 'TOWERINFO_CALIB_HCALOUT'" << std::endl;
0110     return Fun4AllReturnCodes::ABORTRUN;
0111   }
0112 
0113   m_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0114   if (!m_geomOH)
0115   {
0116     std::cerr << Name() << ": could not find RawTowerGeomContainer node 'TOWERGEOM_HCALOUT'" << std::endl;
0117     return Fun4AllReturnCodes::ABORTRUN;
0118   }
0119 
0120   CreateNodes(topNode);
0121   return Fun4AllReturnCodes::EVENT_OK;
0122 }
0123 
0124 void PhotonClusterBuilder::CreateNodes(PHCompositeNode* topNode)
0125 {
0126   PHNodeIterator iter(topNode);
0127   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0128   if (!dstNode)
0129   {
0130     throw std::runtime_error("PhotonClusterBuilder: DST node not found");
0131   }
0132   m_photon_container = findNode::getClass<RawClusterContainer>(dstNode, m_output_photon_node);
0133   if (!m_photon_container)
0134   {
0135     m_photon_container = new RawClusterContainer();
0136     auto* photonNode = new PHIODataNode<PHObject>(m_photon_container, m_output_photon_node, "PHObject");
0137     dstNode->addNode(photonNode);
0138   }
0139 }
0140 
0141 int PhotonClusterBuilder::process_event(PHCompositeNode* topNode)
0142 {
0143   if (!m_rawclusters)
0144   {
0145     m_rawclusters = findNode::getClass<RawClusterContainer>(topNode, m_input_cluster_node);
0146     if (!m_rawclusters)
0147     {
0148       std::cerr << Name() << ": missing RawClusterContainer '" << m_input_cluster_node << "'" << std::endl;
0149       return Fun4AllReturnCodes::ABORTEVENT;
0150     }
0151   }
0152 
0153   // init with NaN
0154   m_vertex = std::numeric_limits<float>::quiet_NaN();
0155   // assume we need vertex for photon shower shape for now
0156   // in the future we need to change the vertex to MBD tracking combined
0157   MbdVertexMap* vertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0158 
0159   if (!vertexmap)
0160   {
0161     std::cout << "GlobalVertexMap node is missing" << std::endl;
0162     return Fun4AllReturnCodes::EVENT_OK;
0163   }
0164   if (vertexmap && !vertexmap->empty())
0165   {
0166     MbdVertex* vtx = vertexmap->begin()->second;
0167     if (vtx)
0168     {
0169       m_vertex = vtx->get_z();
0170 
0171       if (m_vertex != m_vertex)
0172       {
0173         return Fun4AllReturnCodes::EVENT_OK;
0174       }
0175     }
0176     else
0177     {
0178       return Fun4AllReturnCodes::EVENT_OK;
0179     }
0180   }
0181   else
0182   {
0183     return Fun4AllReturnCodes::EVENT_OK;
0184   }
0185 
0186   // iterate over clusters via map to have access to keys if needed
0187   const auto& rcmap = m_rawclusters->getClustersMap();
0188   for (const auto& kv : rcmap)
0189   {
0190     RawCluster* rc = kv.second;
0191     if (!rc)
0192     {
0193       continue;
0194     }
0195 
0196     CLHEP::Hep3Vector vertex_vec(0, 0, m_vertex);
0197 
0198 
0199 
0200     float eta = RawClusterUtility::GetPseudorapidity(*rc, vertex_vec);
0201     float phi = RawClusterUtility::GetAzimuthAngle(*rc, vertex_vec);
0202     float E = rc->get_energy();
0203     float ET = E / std::cosh(eta);
0204     if (ET < m_min_cluster_et)
0205     {
0206       continue;
0207     }
0208 
0209     PhotonClusterv1* photon = new PhotonClusterv1(*rc);
0210 
0211     calculate_shower_shapes(rc, photon, eta, phi);
0212     //this is defensive coding, if do bdt is set false the bdt object should be nullptr
0213     //and this method will simply pass
0214     if (m_do_bdt)
0215     {
0216       calculate_bdt_score(photon);
0217     }
0218 
0219     m_photon_container->AddCluster(photon);
0220   }
0221   return Fun4AllReturnCodes::EVENT_OK;
0222 }
0223 
0224 void PhotonClusterBuilder::calculate_bdt_score(PhotonClusterv1* photon)
0225 {
0226   if (!m_bdt)
0227   {
0228     return;
0229   }
0230 
0231   std::vector<float> x;
0232   for (const auto& feature : m_bdt_feature_list)
0233   {
0234     if (feature == "e11_over_e33")
0235     {
0236       float e11 = photon->get_shower_shape_parameter("e11");
0237       float e33 = photon->get_shower_shape_parameter("e33");
0238       x.push_back((e33 > 0) ? e11 / e33 : 0);
0239     }
0240     else if (feature == "e32_over_e35")
0241     {
0242       float e32 = photon->get_shower_shape_parameter("e32");
0243       float e35 = photon->get_shower_shape_parameter("e35");
0244       x.push_back((e35 > 0) ? e32 / e35 : 0);
0245     }
0246     else if (feature == "vertex_z")
0247     {
0248       x.push_back(m_vertex);
0249     }
0250     else if (feature == "ET")
0251     {
0252       float E = photon->get_energy();
0253       float ET = E / std::cosh(photon->get_shower_shape_parameter("cluster_eta"));
0254       x.push_back(ET);
0255     }
0256     else
0257     {
0258       x.push_back(photon->get_shower_shape_parameter(feature));
0259     }
0260   }
0261 
0262   float bdt_score = -1;  // default value
0263 
0264   bdt_score = m_bdt->Compute(x)[0];
0265 
0266   photon->set_shower_shape_parameter("bdt_score", bdt_score);
0267 }
0268 
0269 void PhotonClusterBuilder::calculate_shower_shapes(RawCluster* rc, PhotonClusterv1* photon, float cluster_eta, float cluster_phi)
0270 {
0271   std::vector<float> showershape = rc->get_shower_shapes(m_shape_min_tower_E);
0272   if (showershape.empty())
0273   {
0274     return;
0275   }
0276 
0277   std::pair<int, int> leadtowerindex = rc->get_lead_tower();
0278   int lead_ieta = leadtowerindex.first;
0279   int lead_iphi = leadtowerindex.second;
0280 
0281   float avg_eta = showershape[4] + 0.5F;
0282   float avg_phi = showershape[5] + 0.5F;
0283 
0284 
0285   int maxieta = std::floor(avg_eta);
0286   int maxiphi = std::floor(avg_phi);
0287 
0288   //if (maxieta < 3 || maxieta > 92)
0289   //{
0290   //  return;
0291   //}
0292 
0293   // for detamax, dphimax, nsaturated
0294   int detamax = 0;
0295   int dphimax = 0;
0296   int nsaturated = 0;
0297   float clusteravgtime = 0;
0298   float cluster_total_e = 0;
0299   const RawCluster::TowerMap& tower_map = rc->get_towermap();
0300   std::set<unsigned int> towers_in_cluster;
0301   for (auto tower_iter : tower_map)
0302   {
0303     RawTowerDefs::keytype tower_key = tower_iter.first;
0304     int ieta = RawTowerDefs::decode_index1(tower_key);
0305     int iphi = RawTowerDefs::decode_index2(tower_key);
0306 
0307     
0308     unsigned int towerinfokey = TowerInfoDefs::encode_emcal(ieta, iphi);
0309     towers_in_cluster.insert(towerinfokey);
0310     TowerInfo* towerinfo = m_emc_tower_container->get_tower_at_key(towerinfokey);
0311     if (towerinfo)
0312     {
0313       clusteravgtime += towerinfo->get_time() * towerinfo->get_energy();
0314       cluster_total_e += towerinfo->get_energy();
0315       if (towerinfo->get_isSaturated())
0316       {
0317         nsaturated++;
0318       }
0319     }
0320     
0321     int totalphibins = 256;
0322     auto dphiwrap = [totalphibins](int towerphi, int maxiphi_arg)
0323     {
0324       int idphi = towerphi - maxiphi_arg;
0325       if (idphi > totalphibins / 2)
0326       {
0327         idphi -= totalphibins;
0328       }
0329       if (idphi < -totalphibins / 2)
0330       {
0331         idphi += totalphibins;
0332       }
0333       return idphi;
0334     };
0335 
0336     int deta = ieta - lead_ieta;
0337     int dphi_val = dphiwrap(iphi, lead_iphi);
0338 
0339     detamax = std::max(std::abs(deta), detamax);
0340     dphimax = std::max(std::abs(dphi_val), dphimax);
0341   }
0342 
0343   if (cluster_total_e > 0)
0344   {
0345     clusteravgtime /= cluster_total_e;
0346   }
0347   else
0348   {
0349     std::cout << "cluster_total_e is 0(this should not happen!!!), setting clusteravgtime to NaN" << std::endl;
0350     clusteravgtime = std::numeric_limits<float>::quiet_NaN();
0351   }
0352 
0353   float E77[7][7] = {{0.0F}};
0354   int E77_ownership[7][7] = {{0}};
0355 
0356   for (int ieta = maxieta - 3; ieta < maxieta + 4; ieta++)
0357   {
0358     for (int iphi = maxiphi - 3; iphi < maxiphi + 4; iphi++)
0359     {
0360       //this is defensive coding, if ieta is out of range, set the energy to 0
0361       //even without this, the requirement for towerinfo object will take care of it
0362       if (ieta < 0 || ieta > 95)
0363       {
0364         E77[ieta - maxieta + 3][iphi - maxiphi + 3] = 0.0F;
0365         E77_ownership[ieta - maxieta + 3][iphi - maxiphi + 3] = 0;
0366         continue;
0367       }
0368 
0369       int temp_ieta = ieta;
0370       int temp_iphi = iphi;
0371       shift_tower_index(temp_ieta, temp_iphi, 96, 256);
0372       if (temp_ieta < 0)
0373       {
0374         continue;
0375       }
0376 
0377       unsigned int towerinfokey = TowerInfoDefs::encode_emcal(temp_ieta, temp_iphi);
0378 
0379       if (towers_in_cluster.contains(towerinfokey))
0380       {
0381         E77_ownership[ieta - maxieta + 3][iphi - maxiphi + 3] = 1;
0382       }
0383 
0384       TowerInfo* towerinfo = m_emc_tower_container->get_tower_at_key(towerinfokey);
0385       if (towerinfo && towerinfo->get_isGood())
0386       {
0387         float energy = towerinfo->get_energy();
0388         if (energy > m_shape_min_tower_E)
0389         {
0390           E77[ieta - maxieta + 3][iphi - maxiphi + 3] = energy;
0391         }
0392       }
0393     }
0394   }
0395 
0396   float e11 = E77[3][3];
0397   float e33 = 0;
0398   float e55 = 0;
0399   float e77 = 0;
0400   float e13 = 0;
0401   float e15 = 0;
0402   float e17 = 0;
0403   float e31 = 0;
0404   float e51 = 0;
0405   float e71 = 0;
0406   float e35 = 0;
0407   float e37 = 0;
0408   float e53 = 0;
0409   float e73 = 0;
0410   float e57 = 0;
0411   float e75 = 0;
0412   float weta = 0;
0413   float wphi = 0;
0414   float weta_cog = 0;
0415   float wphi_cog = 0;
0416   float weta_cogx = 0;
0417   float wphi_cogx = 0;
0418   float Eetaphi = 0;
0419   float shift_eta = avg_eta - std::floor(avg_eta) - 0.5;
0420   float shift_phi = avg_phi - std::floor(avg_phi) - 0.5;
0421   float cog_eta = 3 + shift_eta;
0422   float cog_phi = 3 + shift_phi;
0423 
0424   float w32 = 0;
0425   float e32 = 0;
0426   float w52 = 0;
0427   float e52 = 0;
0428   float w72 = 0;
0429   float e72 = 0;
0430   float detacog = std::abs(maxieta - avg_eta);
0431   float dphicog = std::abs(maxiphi - avg_phi);
0432   float drad = std::sqrt(dphicog*dphicog + detacog*detacog);
0433 
0434 
0435   int signphi = (avg_phi - std::floor(avg_phi)) > 0.5 ? 1 : -1;
0436 
0437   for (int i = 0; i < 7; i++)
0438   {
0439     for (int j = 0; j < 7; j++)
0440     {
0441       int di = std::abs(i - 3);
0442       int dj = std::abs(j - 3);
0443       float di_float = i - cog_eta;
0444       float dj_float = j - cog_phi;
0445 
0446       if (E77_ownership[i][j] == 1)
0447       {
0448         weta += E77[i][j] * di * di;
0449         wphi += E77[i][j] * dj * dj;
0450         weta_cog += E77[i][j] * di_float * di_float;
0451         wphi_cog += E77[i][j] * dj_float * dj_float;
0452         Eetaphi += E77[i][j];
0453         if (i != 3 || j != 3)
0454         {
0455           weta_cogx += E77[i][j] * di_float * di_float;
0456           wphi_cogx += E77[i][j] * dj_float * dj_float;
0457         }
0458       }
0459 
0460       e77 += E77[i][j];
0461       if (di <= 1 && (dj == 0 || j == (3 + signphi)))
0462       {
0463         w32 += E77[i][j] * (i - 3) * (i - 3);
0464         e32 += E77[i][j];
0465       }
0466       if (di <= 2 && (dj == 0 || j == (3 + signphi)))
0467       {
0468         w52 += E77[i][j] * (i - 3) * (i - 3);
0469         e52 += E77[i][j];
0470       }
0471       if (di <= 3 && (dj == 0 || j == (3 + signphi)))
0472       {
0473         w72 += E77[i][j] * (i - 3) * (i - 3);
0474         e72 += E77[i][j];
0475       }
0476 
0477       if (di <= 0 && dj <= 1)
0478       {
0479         e13 += E77[i][j];
0480       }
0481       if (di <= 0 && dj <= 2)
0482       {
0483         e15 += E77[i][j];
0484       }
0485       if (di <= 0 && dj <= 3)
0486       {
0487         e17 += E77[i][j];
0488       }
0489       if (di <= 1 && dj <= 0)
0490       {
0491         e31 += E77[i][j];
0492       }
0493       if (di <= 2 && dj <= 0)
0494       {
0495         e51 += E77[i][j];
0496       }
0497       if (di <= 3 && dj <= 0)
0498       {
0499         e71 += E77[i][j];
0500       }
0501       if (di <= 1 && dj <= 1)
0502       {
0503         e33 += E77[i][j];
0504       }
0505       if (di <= 1 && dj <= 2)
0506       {
0507         e35 += E77[i][j];
0508       }
0509       if (di <= 1 && dj <= 3)
0510       {
0511         e37 += E77[i][j];
0512       }
0513       if (di <= 2 && dj <= 1)
0514       {
0515         e53 += E77[i][j];
0516       }
0517       if (di <= 3 && dj <= 1)
0518       {
0519         e73 += E77[i][j];
0520       }
0521       if (di <= 2 && dj <= 2)
0522       {
0523         e55 += E77[i][j];
0524       }
0525       if (di <= 2 && dj <= 3)
0526       {
0527         e57 += E77[i][j];
0528       }
0529       if (di <= 3 && dj <= 2)
0530       {
0531         e75 += E77[i][j];
0532       }
0533     }
0534   }
0535 
0536   if (Eetaphi > 0)
0537   {
0538     weta /= Eetaphi;
0539     wphi /= Eetaphi;
0540     weta_cog /= Eetaphi;
0541     wphi_cog /= Eetaphi;
0542     weta_cogx /= Eetaphi;
0543     wphi_cogx /= Eetaphi;
0544   }
0545   if (e32 > 0)
0546   {
0547     w32 /= e32;
0548   }
0549   if (e52 > 0)
0550   {
0551     w52 /= e52;
0552   }
0553   if (e72 > 0)
0554   {
0555     w72 /= e72;
0556   }
0557 
0558   photon->set_shower_shape_parameter("et1", showershape[0]);
0559   photon->set_shower_shape_parameter("et2", showershape[1]);
0560   photon->set_shower_shape_parameter("et3", showershape[2]);
0561   photon->set_shower_shape_parameter("et4", showershape[3]);
0562   photon->set_shower_shape_parameter("e11", e11);
0563   photon->set_shower_shape_parameter("e33", e33);
0564   photon->set_shower_shape_parameter("e55", e55);
0565   photon->set_shower_shape_parameter("e77", e77);
0566   photon->set_shower_shape_parameter("e13", e13);
0567   photon->set_shower_shape_parameter("e15", e15);
0568   photon->set_shower_shape_parameter("e17", e17);
0569   photon->set_shower_shape_parameter("e31", e31);
0570   photon->set_shower_shape_parameter("e51", e51);
0571   photon->set_shower_shape_parameter("e71", e71);
0572   photon->set_shower_shape_parameter("e35", e35);
0573   photon->set_shower_shape_parameter("e37", e37);
0574   photon->set_shower_shape_parameter("e53", e53);
0575   photon->set_shower_shape_parameter("e73", e73);
0576   photon->set_shower_shape_parameter("e57", e57);
0577   photon->set_shower_shape_parameter("e75", e75);
0578   photon->set_shower_shape_parameter("weta", weta);
0579   photon->set_shower_shape_parameter("wphi", wphi);
0580   photon->set_shower_shape_parameter("weta_cog", weta_cog);
0581   photon->set_shower_shape_parameter("wphi_cog", wphi_cog);
0582   photon->set_shower_shape_parameter("weta_cogx", weta_cogx);
0583   photon->set_shower_shape_parameter("wphi_cogx", wphi_cogx);
0584   photon->set_shower_shape_parameter("detamax", detamax);
0585   photon->set_shower_shape_parameter("dphimax", dphimax);
0586   photon->set_shower_shape_parameter("nsaturated", nsaturated);
0587   photon->set_shower_shape_parameter("e32", e32);
0588   photon->set_shower_shape_parameter("e52", e52);
0589   photon->set_shower_shape_parameter("e72", e72);
0590   photon->set_shower_shape_parameter("w32", w32);
0591   photon->set_shower_shape_parameter("w52", w52);
0592   photon->set_shower_shape_parameter("w72", w72);
0593   photon->set_shower_shape_parameter("cluster_eta", cluster_eta);
0594   photon->set_shower_shape_parameter("cluster_phi", cluster_phi);
0595   photon->set_shower_shape_parameter("mean_time", clusteravgtime);
0596   photon->set_shower_shape_parameter("detacog", detacog);
0597   photon->set_shower_shape_parameter("dphicog", dphicog);
0598   photon->set_shower_shape_parameter("drad", drad);
0599 
0600   // HCAL info
0601   std::vector<int> ihcal_tower = find_closest_hcal_tower(cluster_eta, cluster_phi, m_geomIH, m_ihcal_tower_container, 0.0, true);
0602   std::vector<int> ohcal_tower = find_closest_hcal_tower(cluster_eta, cluster_phi, m_geomOH, m_ohcal_tower_container, 0.0, false);
0603 
0604   float ihcal_et = 0;
0605   float ohcal_et = 0;
0606   float ihcal_et22 = 0;
0607   float ohcal_et22 = 0;
0608   float ihcal_et33 = 0;
0609   float ohcal_et33 = 0;
0610 
0611   int ihcal_ieta = ihcal_tower[0];
0612   int ihcal_iphi = ihcal_tower[1];
0613   float ihcalEt33[3][3] = {{0.0F}};
0614 
0615   int ohcal_ieta = ohcal_tower[0];
0616   int ohcal_iphi = ohcal_tower[1];
0617   float ohcalEt33[3][3] = {{0.0F}};
0618 
0619   for (int ieta_h = ihcal_ieta - 1; ieta_h <= ihcal_ieta + 1; ieta_h++)
0620   {
0621     for (int iphi_h = ihcal_iphi - 1; iphi_h <= ihcal_iphi + 1; iphi_h++)
0622     {
0623       int temp_ieta = ieta_h;
0624       int temp_iphi = iphi_h;
0625       shift_tower_index(temp_ieta, temp_iphi, 24, 64);
0626       if (temp_ieta < 0)
0627       {
0628         continue;
0629       }
0630 
0631       unsigned int towerinfokey = TowerInfoDefs::encode_hcal(temp_ieta, temp_iphi);
0632       TowerInfo* towerinfo = m_ihcal_tower_container->get_tower_at_key(towerinfokey);
0633       if (towerinfo && towerinfo->get_isGood())
0634       {
0635         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, temp_ieta, temp_iphi);
0636         RawTowerGeom* tower_geom = m_geomIH->get_tower_geometry(key);
0637         if (tower_geom)
0638         {
0639           float energy = towerinfo->get_energy();
0640           float eta = getTowerEta(tower_geom, 0, 0, 0);
0641           float sintheta = 1.0 / std::cosh(eta);
0642           float Et = energy * sintheta;
0643           ihcalEt33[ieta_h - ihcal_ieta + 1][iphi_h - ihcal_iphi + 1] = Et;
0644         }
0645       }
0646     }
0647   }
0648 
0649   for (int ieta_h = ohcal_ieta - 1; ieta_h <= ohcal_ieta + 1; ieta_h++)
0650   {
0651     for (int iphi_h = ohcal_iphi - 1; iphi_h <= ohcal_iphi + 1; iphi_h++)
0652     {
0653       int temp_ieta = ieta_h;
0654       int temp_iphi = iphi_h;
0655       shift_tower_index(temp_ieta, temp_iphi, 24, 64);
0656       if (temp_ieta < 0)
0657       {
0658         continue;
0659       }
0660 
0661       unsigned int towerinfokey = TowerInfoDefs::encode_hcal(temp_ieta, temp_iphi);
0662       TowerInfo* towerinfo = m_ohcal_tower_container->get_tower_at_key(towerinfokey);
0663       if (towerinfo && towerinfo->get_isGood())
0664       {
0665         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, temp_ieta, temp_iphi);
0666         RawTowerGeom* tower_geom = m_geomOH->get_tower_geometry(key);
0667         if (tower_geom)
0668         {
0669           float energy = towerinfo->get_energy();
0670           float eta = getTowerEta(tower_geom, 0, 0, 0);
0671           float sintheta = 1.0 / std::cosh(eta);
0672           float Et = energy * sintheta;
0673           ohcalEt33[ieta_h - ohcal_ieta + 1][iphi_h - ohcal_iphi + 1] = Et;
0674         }
0675       }
0676     }
0677   }
0678 
0679   ihcal_et = ihcalEt33[1][1];
0680   ohcal_et = ohcalEt33[1][1];
0681 
0682   for (int i = 0; i < 3; i++)
0683   {
0684     for (int j = 0; j < 3; j++)
0685     {
0686       ihcal_et33 += ihcalEt33[i][j];
0687       ohcal_et33 += ohcalEt33[i][j];
0688       if (i == 1 || j == 1 + ihcal_tower[2])
0689       {
0690         if (j == 1 || i == 1 + ihcal_tower[3])
0691         {
0692           ihcal_et22 += ihcalEt33[i][j];
0693         }
0694       }
0695       if (i == 1 || j == 1 + ohcal_tower[2])
0696       {
0697         if (j == 1 || i == 1 + ohcal_tower[3])
0698         {
0699           ohcal_et22 += ohcalEt33[i][j];
0700         }
0701       }
0702     }
0703   }
0704 
0705   photon->set_shower_shape_parameter("ihcal_et", ihcal_et);
0706   photon->set_shower_shape_parameter("ohcal_et", ohcal_et);
0707   photon->set_shower_shape_parameter("ihcal_et22", ihcal_et22);
0708   photon->set_shower_shape_parameter("ohcal_et22", ohcal_et22);
0709   photon->set_shower_shape_parameter("ihcal_et33", ihcal_et33);
0710   photon->set_shower_shape_parameter("ohcal_et33", ohcal_et33);
0711   photon->set_shower_shape_parameter("ihcal_ieta", ihcal_ieta);
0712   photon->set_shower_shape_parameter("ihcal_iphi", ihcal_iphi);
0713   photon->set_shower_shape_parameter("ohcal_ieta", ohcal_ieta);
0714   photon->set_shower_shape_parameter("ohcal_iphi", ohcal_iphi);
0715 
0716   float E = photon->get_energy();
0717   float ET = E / std::cosh(cluster_eta);
0718 
0719   auto compute_layer_iso = [&](RawTowerDefs::CalorimeterId calo_id, float radius)
0720   {
0721     TowerInfoContainer* container = nullptr;
0722     RawTowerGeomContainer* geom = nullptr;
0723     if (calo_id == RawTowerDefs::CalorimeterId::CEMC)
0724     {
0725       container = m_emc_tower_container;
0726       geom = m_geomEM;
0727     }
0728     else if (calo_id == RawTowerDefs::CalorimeterId::HCALIN)
0729     {
0730       container = m_ihcal_tower_container;
0731       geom = m_geomIH;
0732     }
0733     else
0734     {
0735       container = m_ohcal_tower_container;
0736       geom = m_geomOH;
0737     }
0738     return calculate_layer_et(cluster_eta, cluster_phi, radius, container, geom, calo_id, m_vertex);
0739   };
0740 
0741   const float emcal_et_04 = compute_layer_iso(RawTowerDefs::CalorimeterId::CEMC, 0.4);
0742   const float ihcal_et_04 = compute_layer_iso(RawTowerDefs::CalorimeterId::HCALIN, 0.4);
0743   const float ohcal_et_04 = compute_layer_iso(RawTowerDefs::CalorimeterId::HCALOUT, 0.4);
0744 
0745   const float emcal_et_03 = compute_layer_iso(RawTowerDefs::CalorimeterId::CEMC, 0.3);
0746   const float ihcal_et_03 = compute_layer_iso(RawTowerDefs::CalorimeterId::HCALIN, 0.3);
0747   const float ohcal_et_03 = compute_layer_iso(RawTowerDefs::CalorimeterId::HCALOUT, 0.3);
0748 
0749   const float emcal_et_02 = compute_layer_iso(RawTowerDefs::CalorimeterId::CEMC, 0.2);
0750   const float emcal_et_01 = compute_layer_iso(RawTowerDefs::CalorimeterId::CEMC, 0.1);
0751   const float emcal_et_005 = compute_layer_iso(RawTowerDefs::CalorimeterId::CEMC, 0.05);
0752 
0753   photon->set_shower_shape_parameter("iso_04_emcal", emcal_et_04 - ET);
0754   photon->set_shower_shape_parameter("iso_04_hcalin", ihcal_et_04);
0755   photon->set_shower_shape_parameter("iso_04_hcalout", ohcal_et_04);
0756 
0757   photon->set_shower_shape_parameter("iso_03_emcal", emcal_et_03 - ET);
0758   photon->set_shower_shape_parameter("iso_03_hcalin", ihcal_et_03);
0759   photon->set_shower_shape_parameter("iso_03_hcalout", ohcal_et_03);
0760   photon->set_shower_shape_parameter("iso_02_emcal", emcal_et_02 - ET);
0761   photon->set_shower_shape_parameter("iso_01_emcal", emcal_et_01 - ET);
0762   photon->set_shower_shape_parameter("iso_005_emcal", emcal_et_005 - ET);
0763 }
0764 
0765 double PhotonClusterBuilder::getTowerEta(RawTowerGeom* tower_geom, double vx, double vy, double vz)
0766 {
0767   if (!tower_geom)
0768   {
0769     return -9999;
0770   }
0771   if (vx == 0 && vy == 0 && vz == 0)
0772   {
0773     return tower_geom->get_eta();
0774   }
0775 
0776   double radius = sqrt((tower_geom->get_center_x() - vx) * (tower_geom->get_center_x() - vx) + (tower_geom->get_center_y() - vy) * (tower_geom->get_center_y() - vy));
0777   double theta = atan2(radius, tower_geom->get_center_z() - vz);
0778   return -log(tan(theta / 2.));
0779 }
0780 
0781 std::vector<int> PhotonClusterBuilder::find_closest_hcal_tower(float eta, float phi, RawTowerGeomContainer* geom, TowerInfoContainer* towerContainer, float vertex_z, bool isihcal)
0782 {
0783   int matchedieta = -1;
0784   int matchediphi = -1;
0785   double matchedeta = -999;
0786   double matchedphi = -999;
0787 
0788   if (!geom || !towerContainer)
0789   {
0790     return {-1, -1, 0, 0};
0791   }
0792 
0793   unsigned int ntowers = towerContainer->size();
0794   float minR = 999;
0795 
0796   for (unsigned int channel = 0; channel < ntowers; channel++)
0797   {
0798     TowerInfo* tower = towerContainer->get_tower_at_channel(channel);
0799     if (!tower)
0800     {
0801       continue;
0802     }
0803 
0804     unsigned int towerkey = towerContainer->encode_key(channel);
0805     int ieta = towerContainer->getTowerEtaBin(towerkey);
0806     int iphi = towerContainer->getTowerPhiBin(towerkey);
0807 
0808     RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(isihcal ? RawTowerDefs::CalorimeterId::HCALIN : RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0809     RawTowerGeom* tower_geom = geom->get_tower_geometry(key);
0810     if (!tower_geom)
0811     {
0812       continue;
0813     }
0814 
0815     double this_phi = tower_geom->get_phi();
0816     double this_eta = getTowerEta(tower_geom, 0, 0, vertex_z);
0817     double dR_val = deltaR(eta, phi, this_eta, this_phi);
0818     if (dR_val < minR)
0819     {
0820       minR = dR_val;
0821       matchedieta = ieta;
0822       matchediphi = iphi;
0823       matchedeta = this_eta;
0824       matchedphi = this_phi;
0825     }
0826   }
0827 
0828   float deta = eta - matchedeta;
0829   float dphi_val = phi - matchedphi;
0830   if (dphi_val > M_PI)
0831   {
0832     dphi_val -= 2 * M_PI;
0833   }
0834   if (dphi_val < -M_PI)
0835   {
0836     dphi_val += 2 * M_PI;
0837   }
0838 
0839   int dphisign = (dphi_val > 0) ? 1 : -1;
0840   int detasign = (deta > 0) ? 1 : -1;
0841 
0842   return {matchedieta, matchediphi, detasign, dphisign};
0843 }
0844 
0845 float PhotonClusterBuilder::calculate_layer_et(float seed_eta, float seed_phi, float radius, TowerInfoContainer* towerContainer, RawTowerGeomContainer* geomContainer, RawTowerDefs::CalorimeterId calo_id, float vertex_z)
0846 {
0847   if (!towerContainer || !geomContainer)
0848   {
0849     return std::numeric_limits<float>::quiet_NaN();
0850   }
0851 
0852   float layer_et = 0.0;
0853   const unsigned int ntowers = towerContainer->size();
0854   for (unsigned int channel = 0; channel < ntowers; ++channel)
0855   {
0856     TowerInfo* tower = towerContainer->get_tower_at_channel(channel);
0857     if (!tower || !tower->get_isGood())
0858     {
0859       continue;
0860     }
0861 
0862     unsigned int towerkey = towerContainer->encode_key(channel);
0863     int ieta = towerContainer->getTowerEtaBin(towerkey);
0864     int iphi = towerContainer->getTowerPhiBin(towerkey);
0865 
0866     RawTowerDefs::keytype geom_key = RawTowerDefs::encode_towerid(calo_id, ieta, iphi);
0867     RawTowerGeom* tower_geom = geomContainer->get_tower_geometry(geom_key);
0868     if (!tower_geom)
0869     {
0870       continue;
0871     }
0872 
0873     double tower_eta = getTowerEta(tower_geom, 0, 0, vertex_z);
0874     double tower_phi = tower_geom->get_phi();
0875 
0876     if (deltaR(seed_eta, seed_phi, tower_eta, tower_phi) >= radius)
0877     {
0878       continue;
0879     }
0880 
0881     float energy = tower->get_energy();
0882     if (energy <= m_shape_min_tower_E)
0883     {
0884       continue;
0885     }
0886 
0887     float et = energy / std::cosh(tower_eta);
0888     layer_et += et;
0889   }
0890 
0891   return layer_et;
0892 }
0893 
0894 double PhotonClusterBuilder::deltaR(double eta1, double phi1, double eta2, double phi2)
0895 {
0896   double dphi = phi1 - phi2;
0897   while (dphi > M_PI)
0898   {
0899     dphi -= 2 * M_PI;
0900   }
0901   while (dphi <= -M_PI)
0902   {
0903     dphi += 2 * M_PI;
0904   }
0905   return sqrt(pow(eta1 - eta2, 2) + pow(dphi, 2));
0906 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!