packages/calovtxreco/CaloVtxReco.cc

0001 #include "CaloVtxReco.h"
0002
0003 #include <jetbase/Jet.h>
0004 #include <jetbase/JetContainer.h>
0005
0006 #include <calobase/RawTowerGeom.h>
0007 #include <calobase/RawTowerGeomContainer.h>
0008 #include <calobase/TowerInfo.h>
0009 #include <calobase/TowerInfoContainer.h>
0010 #include <calobase/TowerInfoDefs.h>
0011
0012 #include <globalvertex/CaloVertexMapv1.h>
0013 #include <globalvertex/CaloVertexv1.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/getClass.h>
0019
0020 #include <cmath>
0021 /*
0022 float radius_EM = 93.5;
0023 float radius_OH = 225.87;
0024 */
0025 //____________________________________________________________________________..
0026 CaloVtxReco::CaloVtxReco(const std::string &name, const std::string &jetnodename, const bool use_z_energy_dep)
0027   : SubsysReco(name)
0028   , m_use_z_energy_dep(use_z_energy_dep)
0029   , m_jetnodename(jetnodename)
0030 {
0031 }
0032
0033 int CaloVtxReco::createNodes(PHCompositeNode *topNode)
0034 {
0035   PHNodeIterator iter(topNode);
0036   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0037   if (!dstNode)
0038   {
0039     std::cout << PHWHERE << "DST Node missing doing nothing" << std::endl;
0040     return Fun4AllReturnCodes::ABORTRUN;
0041   }
0042
0043   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0044   if (!runNode)
0045   {
0046     std::cout << PHWHERE << "RUN Node missing doing nothing" << std::endl;
0047     return Fun4AllReturnCodes::ABORTRUN;
0048   }
0049
0050   PHNodeIterator dstiter(dstNode);
0051
0052   PHCompositeNode *globalNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "GLOBAL"));
0053   if (!globalNode)
0054   {
0055     globalNode = new PHCompositeNode("GLOBAL");
0056     dstNode->addNode(globalNode);
0057   }
0058
0059   m_calovtxmap = findNode::getClass<CaloVertexMap>(globalNode, "CaloVertexMap");
0060   if (!m_calovtxmap)
0061   {
0062     m_calovtxmap = new CaloVertexMapv1();
0063     PHIODataNode<PHObject> *VertexMapNode = new PHIODataNode<PHObject>(m_calovtxmap, "CaloVertexMap", "PHObject");
0064     globalNode->addNode(VertexMapNode);
0065   }
0066
0067   return Fun4AllReturnCodes::EVENT_OK;
0068 }
0069
0070 //____________________________________________________________________________..
0071 int CaloVtxReco::InitRun(PHCompositeNode *topNode)
0072 {
0073   if (Verbosity() > 1)
0074   {
0075     std::cout << "Initializing!" << std::endl;
0076   }
0077   if (createNodes(topNode) == Fun4AllReturnCodes::ABORTRUN)
0078   {
0079     return Fun4AllReturnCodes::ABORTRUN;
0080   }
0081
0082   RawTowerGeomContainer *tower_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0083   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0084   if(!tower_geomEM || !tower_geomOH)
0085     {
0086       if(Verbosity() > 0)
0087     {
0088       std::cout << "CaloVtxReco::InitRun(): Missing tower geometry node for towergeomEM (address: " << tower_geomEM << ") or towergeomOH (address: " << tower_geomOH << ") - aborting run!" << std::endl;
0089     }
0090       return Fun4AllReturnCodes::ABORTRUN;
0091     }
0092
0093   m_radius_EM = tower_geomEM->get_radius();
0094   m_radius_OH = tower_geomOH->get_radius();
0095
0096   if(std::isnan(m_radius_EM) || std::isnan(m_radius_OH))
0097     {
0098       if(Verbosity() > 0)
0099     {
0100       std::cout << "CaloVtxReco::InitRun(): NaN value for one of radius EM (value: " << m_radius_EM << ") or radius OH (value: " << m_radius_OH << ") after attempting to get - aborting run!" << std::endl;
0101     }
0102       return Fun4AllReturnCodes::ABORTRUN;
0103     }
0104
0105
0106   return Fun4AllReturnCodes::EVENT_OK;
0107 }
0108
0109 float CaloVtxReco::new_eta(int channel, TowerInfoContainer *towers, RawTowerGeomContainer *geom, RawTowerDefs::CalorimeterId caloID, float testz)
0110 {
0111   int key = towers->encode_key(channel);
0112   const RawTowerDefs::keytype geomkey = RawTowerDefs::encode_towerid(caloID, towers->getTowerEtaBin(key), towers->getTowerPhiBin(key));
0113
0114   RawTowerGeom *tower_geom = geom->get_tower_geometry(geomkey);
0115   float oldeta = tower_geom->get_eta();
0116
0117   float radius = (caloID == RawTowerDefs::CalorimeterId::HCALIN ? m_radius_EM : m_radius_OH);
0118   float towerz = radius / (tanf(2 * std::atanf(std::exp(oldeta))));
0119   float newz = towerz + testz;
0120   float newTheta = std::atan2(radius, newz);
0121   float neweta = -log(tan(0.5 * newTheta));
0122
0123   return neweta;
0124 }
0125
0126 float get_dphi(float phi1, float phi2)
0127 {
0128   float dphi = std::abs(phi1 - phi2);
0129   if (dphi > M_PI)
0130   {
0131     dphi = 2 * M_PI - dphi;
0132   }
0133   return dphi;
0134 }
0135
0136 int CaloVtxReco::calo_tower_algorithm(PHCompositeNode *topNode) const
0137 {
0138   TowerInfoContainer *emcal_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0139   TowerInfoContainer *hcalin_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0140   TowerInfoContainer *hcalout_towers = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0141   RawTowerGeomContainer *tower_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0142   RawTowerGeomContainer *tower_geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0143   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0144
0145   int size;
0146
0147   float average_z[3]{0};
0148   float total_E[3]{0};
0149
0150   if (emcal_towers)
0151   {
0152     size = emcal_towers->size();  // online towers should be the same!
0153     for (int channel = 0; channel < size; channel++)
0154     {
0155       TowerInfo *_tower = emcal_towers->get_tower_at_channel(channel);
0156       short good = (_tower->get_isGood() ? 1 : 0);
0157       if (!good)
0158       {
0159         continue;
0160       }
0161
0162       float energy = _tower->get_energy();
0163       if (energy < m_energy_cut)
0164       {
0165         continue;
0166       }
0167
0168       // float time = _tower->get_time_float();
0169
0170       unsigned int towerkey = emcal_towers->encode_key(channel);
0171       int ieta = emcal_towers->getTowerEtaBin(towerkey);
0172       int iphi = emcal_towers->getTowerPhiBin(towerkey);
0173
0174       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0175       /*
0176       if (emcal_r < 10)
0177         {
0178            emcal_r = tower_geomEM->get_tower_geometry(key)->get_center_radius();
0179         }
0180       */
0181       float tower_z = tower_geomEM->get_tower_geometry(key)->get_center_z();
0182       average_z[0] += tower_z * energy;
0183       total_E[0] += energy;
0184     }
0185   }
0186
0187   if (hcalin_towers)
0188   {
0189     size = hcalin_towers->size();  // online towers should be the same!
0190     for (int channel = 0; channel < size; channel++)
0191     {
0192       TowerInfo *_tower = hcalin_towers->get_tower_at_channel(channel);
0193       float energy = _tower->get_energy();
0194       if (energy < m_energy_cut)
0195       {
0196         continue;
0197       }
0198       // float time = _tower->get_time_float();
0199       short good = (_tower->get_isGood() ? 1 : 0);
0200       if (!good)
0201       {
0202         continue;
0203       }
0204
0205       unsigned int towerkey = hcalin_towers->encode_key(channel);
0206       int ieta = hcalin_towers->getTowerEtaBin(towerkey);
0207       int iphi = hcalin_towers->getTowerPhiBin(towerkey);
0208       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0209       float tower_z = tower_geomIH->get_tower_geometry(key)->get_center_z();
0210       /*
0211       if (hcalin_r < 10)
0212         {
0213            hcalin_r = tower_geomIH->get_tower_geometry(key)->get_center_radius();
0214         }
0215       */
0216       average_z[1] += tower_z * energy;
0217       total_E[1] += energy;
0218     }
0219   }
0220   if (hcalout_towers)
0221   {
0222     size = hcalout_towers->size();  // online towers should be the same!
0223     for (int channel = 0; channel < size; channel++)
0224     {
0225       TowerInfo *_tower = hcalout_towers->get_tower_at_channel(channel);
0226       float energy = _tower->get_energy();
0227       if (energy < m_energy_cut)
0228       {
0229         continue;
0230       }
0231       // float time = _tower->get_time_float();
0232       unsigned int towerkey = hcalout_towers->encode_key(channel);
0233       int ieta = hcalout_towers->getTowerEtaBin(towerkey);
0234       int iphi = hcalout_towers->getTowerPhiBin(towerkey);
0235       short good = (_tower->get_isGood() ? 1 : 0);
0236
0237       if (!good)
0238       {
0239         continue;
0240       }
0241
0242       const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0243       /*
0244       if (hcalout_r < 10)
0245         {
0246           hcalout_r = tower_geomOH->get_tower_geometry(key)->get_center_radius();
0247         }
0248       */
0249       float tower_z = tower_geomOH->get_tower_geometry(key)->get_center_z();
0250       average_z[2] += tower_z * energy;
0251       total_E[2] += energy;
0252     }
0253   }
0254
0255   double b_calo_vertex_z = (average_z[0] + average_z[1] + average_z[2]) / (total_E[0] + total_E[1] + total_E[2]);
0256
0257   CaloVertex *vertex = new CaloVertexv1();
0258   vertex->set_z(b_calo_vertex_z);
0259   m_calovtxmap->insert(vertex);
0260
0261   return 0;
0262 }
0263
0264 int CaloVtxReco::process_event(PHCompositeNode *topNode)
0265 {
0266   if (m_use_z_energy_dep)
0267   {
0268     calo_tower_algorithm(topNode);
0269     return Fun4AllReturnCodes::EVENT_OK;
0270   }
0271
0272   if (Verbosity() > 1)
0273   {
0274     std::cout << std::endl
0275               << std::endl
0276               << std::endl
0277               << "CaloVtxReco: Beginning event processing" << std::endl;
0278   }
0279   m_zvtx = std::numeric_limits<float>::quiet_NaN();
0280   JetContainer *jetcon = findNode::getClass<JetContainer>(topNode, m_jetnodename);
0281   TowerInfoContainer *towers[3];
0282   towers[0] = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0283   towers[1] = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0284   towers[2] = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0285
0286   RawTowerGeomContainer *geom[3];
0287   geom[0] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0288   geom[1] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0289   geom[2] = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0290
0291   const int nz = 601;
0292   const int njet = 2;
0293   Jet *jets[njet];
0294   float jpt[njet] = {0};
0295   float jemsum[njet] = {0};
0296   float johsum[njet] = {0};
0297   float jemeta[njet] = {0};
0298   float joheta[njet] = {0};
0299
0300   if (jetcon)
0301   {
0302     int tocheck = jetcon->size();
0303     if (Verbosity() > 2)
0304     {
0305       std::cout << "Found " << tocheck << " jets to check..." << std::endl;
0306     }
0307     for (int i = 0; i < tocheck; ++i)
0308     {
0309       Jet *jet = jetcon->get_jet(i);
0310       if (jet)
0311       {
0312         float pt = jet->get_pt();
0313         if (pt < m_jet_threshold)
0314         {
0315           continue;
0316         }
0317         if (pt > jpt[0])
0318         {
0319           jpt[1] = jpt[0];
0320           jets[1] = jets[0];
0321           jets[0] = jet;
0322           jpt[0] = pt;
0323         }
0324         else if (pt > jpt[1])
0325         {
0326           jets[1] = jet;
0327           jpt[1] = pt;
0328         }
0329       }
0330     }
0331   }
0332   else
0333   {
0334     if (Verbosity() > 0)
0335     {
0336       std::cout << "no jets" << std::endl;
0337     }
0338     return Fun4AllReturnCodes::ABORTEVENT;
0339   }
0340
0341   if (jpt[0] == 0)
0342   {
0343     if (Verbosity() > 2)
0344     {
0345       std::cout << "NO JETS > 5 GeV!" << std::endl;
0346     }
0347   }
0348
0349   float metric = std::numeric_limits<float>::max();
0350   for (int i = 0; i < nz; ++i)
0351   {
0352     float testz = -300 + i;
0353     float testmetric = 0;
0354     for (int j = 0; j < njet; ++j)
0355     {
0356       if (jpt[j] == 0)
0357       {
0358         continue;
0359       }
0360       jemsum[j] = 0;
0361       johsum[j] = 0;
0362       jemeta[j] = 0;
0363       joheta[j] = 0;
0364       for (auto comp : jets[j]->get_comp_vec())
0365       {
0366         if (comp.first == 5 || comp.first == 26)
0367         {
0368           continue;
0369         }
0370         unsigned int channel = comp.second;
0371         if (comp.first == 7 || comp.first == 27)
0372         {
0373           TowerInfo *tower = towers[2]->get_tower_at_channel(channel);
0374           if (tower->get_energy() < 0.1)
0375           {
0376             continue;
0377           }
0378           johsum[j] += tower->get_energy();
0379           float neweta = new_eta(channel, towers[2], geom[2], RawTowerDefs::CalorimeterId::HCALOUT, testz);
0380           joheta[j] += neweta * tower->get_energy();
0381         }
0382         if (comp.first == 13 || comp.first == 28 || comp.first == 25)
0383         {
0384           TowerInfo *tower = towers[0]->get_tower_at_channel(channel);
0385           if (tower->get_energy() < 0.1)
0386           {
0387             continue;
0388           }
0389           jemsum[j] += tower->get_energy();
0390           float neweta = new_eta(channel, towers[0], geom[1], RawTowerDefs::CalorimeterId::HCALIN, testz);
0391           jemeta[j] += neweta * tower->get_energy();
0392         }
0393       }
0394       jemeta[j] /= jemsum[j];
0395       joheta[j] /= johsum[j];
0396       if ((jemsum[j] == 0 || johsum[j] == 0) && Verbosity() > 1)
0397       {
0398         std::cout << "zero E sum in at least one calo for a jet" << std::endl;
0399       }
0400       if (!std::isnan(jemeta[j]) && !std::isnan(joheta[j]))
0401       {
0402         testmetric += pow(jemeta[j] - joheta[j], 2);
0403       }
0404     }
0405     if (Verbosity() > 3)
0406     {
0407       std::cout << "metric: " << testmetric << std::endl;
0408     }
0409     if (testmetric < metric && testmetric != 0)
0410     {
0411       metric = testmetric;
0412       m_zvtx = testz;
0413     }
0414   }
0415   if (std::abs(m_zvtx) < 305)
0416   {
0417     if (Verbosity() > 2)
0418     {
0419       std::cout << "optimal z: " << m_zvtx << std::endl;
0420     }
0421     CaloVertex *vertex = new CaloVertexv1();
0422     m_zvtx *= m_calib_factor;  // calibration factor from simulation
0423     vertex->set_z(m_zvtx);
0424     m_calovtxmap->insert(vertex);
0425     if (Verbosity() > 3)
0426     {
0427       std::cout << "CaloVtxReco: end event" << std::endl;
0428     }
0429   }
0430   return Fun4AllReturnCodes::EVENT_OK;
0431 }