packages/CaloReco/RawClusterLikelihoodProfile.cc

0001 #include "RawClusterLikelihoodProfile.h"
0002
0003 #include <calobase/RawCluster.h>
0004 #include <calobase/RawClusterContainer.h>
0005 #include <calobase/RawClusterUtility.h>
0006 #include <calobase/RawTower.h>
0007 #include <calobase/RawTowerContainer.h>
0008 #include <calobase/RawTowerDefs.h>
0009 #include <calobase/RawTowerGeom.h>
0010 #include <calobase/RawTowerGeomContainer.h>
0011
0012 #include <calobase/TowerInfo.h>
0013 #include <calobase/TowerInfoContainer.h>
0014 #include <calobase/TowerInfoDefs.h>
0015
0016 #include <globalvertex/GlobalVertex.h>
0017 #include <globalvertex/GlobalVertexMap.h>
0018
0019 #include <ffamodules/CDBInterface.h>
0020 #include <ffaobjects/EventHeader.h>
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/Fun4AllServer.h>
0023 #include <phool/getClass.h>
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/phool.h>
0026
0027 #include <iostream>
0028 #include <vector>
0029
0030 RawClusterLikelihoodProfile::RawClusterLikelihoodProfile(const std::string &name)
0031   : SubsysReco(name)
0032 {
0033 }
0034
0035 int RawClusterLikelihoodProfile::Init(PHCompositeNode *topNode)
0036 {
0037   cdfcalc = new ClusterCDFCalculator();
0038   cdfcalc->LoadProfile(m_profile_path);
0039
0040   if (m_inputNodeName == m_outputNodeName)
0041   {
0042     std::cout << "RawClusterLikelihoodProfile::Init: Same inputNodeName and outputNodeName, putting inplace to true for overwriting" << std::endl;
0043     inplace = true;
0044   }
0045   CreateNodes(topNode);
0046
0047   return Fun4AllReturnCodes::EVENT_OK;
0048 }
0049
0050 int RawClusterLikelihoodProfile::process_event(PHCompositeNode *topNode)
0051 {
0052   std::string clusterNodeName = m_inputNodeName;
0053   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, clusterNodeName);
0054   if (!clusterContainer)
0055   {
0056     std::cout << "RawClusterLikelihoodProfile::process_event::Could not locate input cluster node " << clusterNodeName << std::endl;
0057     return Fun4AllReturnCodes::ABORTEVENT;
0058   }
0059
0060   if (inplace)
0061   {
0062     _clusters = clusterContainer;
0063   }
0064   else
0065   {
0066     _clusters->Reset();
0067     RawClusterContainer::Map clusterMap = clusterContainer->getClustersMap();
0068     for (auto &clusterPair : clusterMap)
0069     {
0070       RawCluster *cluster = (clusterPair.second);
0071       float clusterE = cluster->get_energy();
0072       if (clusterE < m_min_cluster_e)
0073       {
0074         continue;
0075       }
0076       RawCluster *profileCluster = (RawCluster *) cluster->CloneMe();
0077       _clusters->AddCluster(profileCluster);
0078     }
0079   }
0080
0081   std::string towerNodeName = m_towerNodeName;
0082   TowerInfoContainer *emcTowerContainer = findNode::getClass<TowerInfoContainer>(topNode, towerNodeName);
0083   if (!emcTowerContainer)
0084   {
0085     std::cout << "RawClusterLikelihoodProfile::process_event Could not locate tower node " << towerNodeName << std::endl;
0086     return Fun4AllReturnCodes::ABORTEVENT;
0087   }
0088
0089   std::string globalvtxNodeName = m_globalvtxNodeName;
0090   GlobalVertexMap *globalvtxmap = findNode::getClass<GlobalVertexMap>(topNode,globalvtxNodeName);
0091   bool isglbvtx=true;
0092   if(!globalvtxmap || globalvtxmap->empty())
0093   {
0094     isglbvtx=false;
0095   }
0096
0097   if(isglbvtx){
0098     GlobalVertex *bvertex= nullptr;
0099     for (GlobalVertexMap::ConstIter globaliter= globalvtxmap->begin(); globaliter != globalvtxmap->end(); ++globaliter)
0100     {
0101       bvertex = globaliter->second;
0102     }
0103     if(!bvertex){std::cout << "WARNING! RawClusterLikelihoodProfile::could not find globalvtxmap iter :: set vtx to (0,0,0)" << std::endl;}
0104     else if(bvertex){
0105       vz = bvertex->get_z();
0106       vy = bvertex->get_y();
0107       vx = bvertex->get_x();
0108     }
0109   }
0110
0111   RawClusterContainer::Map clusterMap = _clusters->getClustersMap();
0112   for (auto &clusterPair : clusterMap)
0113   {
0114     RawCluster *cluster = clusterPair.second;
0115     cluster->set_prob(-1);
0116     cluster->set_merged_cluster_prob(-1);
0117     CLHEP::Hep3Vector vertex(vx, vy, vz);
0118     CLHEP::Hep3Vector E_vec_cluster_Full = RawClusterUtility::GetEVec(*cluster, vertex);
0119     double clusterE = E_vec_cluster_Full.mag();
0120     if (clusterE < m_min_cluster_e)
0121     {
0122       continue;
0123     }
0124     const RawCluster::TowerMap tower_map = cluster->get_towermap();
0125
0126     std::vector<float> shower_shapes = cluster->get_shower_shapes(m_tower_thres_e);
0127     int ieta_center_of_gravity = std::floor(shower_shapes[4] + 0.5);
0128     int iphi_center_of_gravity = std::floor(shower_shapes[5] + 0.5);
0129
0130     std::vector<double> input;
0131     int vectorSize = inputDimx * inputDimy;
0132     input.resize(vectorSize, 0);
0133
0134     int xlength = ((inputDimx - 1) / 2);
0135     int ylength = ((inputDimy - 1) / 2);
0136     if (ieta_center_of_gravity - ylength < 0 || ieta_center_of_gravity + ylength >= 96)
0137     {
0138       continue;
0139     }
0140     for (int ieta = ieta_center_of_gravity - ylength; ieta <= ieta_center_of_gravity + ylength; ieta++)
0141     {
0142       for (int iphi = iphi_center_of_gravity - xlength; iphi <= iphi_center_of_gravity + xlength; iphi++)
0143       {
0144         int mappediphi = iphi;
0145
0146         if (mappediphi < 0)
0147         {
0148           mappediphi += 256;
0149         }
0150         if (mappediphi > 255)
0151         {
0152           mappediphi -= 256;
0153         }
0154         unsigned int towerinfokey = TowerInfoDefs::encode_emcal(ieta, mappediphi);
0155         TowerInfo *towerinfo = emcTowerContainer->get_tower_at_key(towerinfokey);
0156         if (!towerinfo)
0157         {
0158           std::cout << "No towerinfo for tower key " << towerinfokey << std::endl;
0159           std::cout << "ieta: " << ieta << " iphi: " << mappediphi << std::endl;
0160           continue;
0161         }
0162         int index = ((ieta - ieta_center_of_gravity + ylength) * inputDimx) + iphi - iphi_center_of_gravity + xlength;
0163         input.at(index) = towerinfo->get_energy();
0164       }
0165     }
0166     std::pair<double, double> probpair = cdfcalc->GetCDF(input, clusterE, m_profile_dimension);
0167     double prob = probpair.first;
0168     double prob_merged_cluster = probpair.second;
0169     cluster->set_prob(prob);
0170     cluster->set_merged_cluster_prob(prob_merged_cluster);
0171   }
0172
0173   return Fun4AllReturnCodes::EVENT_OK;
0174 }
0175
0176 int RawClusterLikelihoodProfile::End(PHCompositeNode * /*topNode*/)
0177 {
0178   delete cdfcalc;
0179   return Fun4AllReturnCodes::EVENT_OK;
0180 }
0181
0182 void RawClusterLikelihoodProfile::Print(const std::string &what) const
0183 {
0184   std::cout << "RawClusterLikelihoodProfile::Print(const std::string &what) const Printing info for " << what << std::endl;
0185   return;
0186 }
0187
0188 void RawClusterLikelihoodProfile::CreateNodes(PHCompositeNode *topNode)
0189 {
0190   PHNodeIterator iter(topNode);
0191
0192   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0193   if (!dstNode)
0194   {
0195     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0196     throw std::runtime_error("Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0197   }
0198
0199   PHCompositeNode *cemcNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "CEMC"));
0200
0201   if (!cemcNode)
0202   {
0203     cemcNode = new PHCompositeNode("CEMC");
0204     dstNode->addNode(cemcNode);
0205   }
0206   std::string clusterNodeName = m_outputNodeName;
0207   _clusters = findNode::getClass<RawClusterContainer>(dstNode, clusterNodeName);
0208   if (!_clusters)
0209   {
0210     _clusters = new RawClusterContainer();
0211     PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_clusters, clusterNodeName, "PHObject");
0212     cemcNode->addNode(clusterNode);
0213   }
0214 }