g4simulation/g4calo/RawTowerBuilderByHitIndex.cc

0001 #include "RawTowerBuilderByHitIndex.h"
0002
0003 #include <calobase/RawTowerContainer.h>
0004 #include <calobase/RawTowerv1.h>
0005
0006 #include <calobase/RawTower.h>               // for RawTower
0007 #include <calobase/RawTowerDefs.h>           // for convert_name_to_caloid
0008 #include <calobase/RawTowerGeom.h>           // for RawTowerGeom
0009 #include <calobase/RawTowerGeomContainer.h>  // for RawTowerGeomContainer
0010 #include <calobase/RawTowerGeomContainerv1.h>
0011 #include <calobase/RawTowerGeomv3.h>
0012
0013 #include <g4main/PHG4Hit.h>
0014 #include <g4main/PHG4HitContainer.h>
0015
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017 #include <fun4all/SubsysReco.h>  // for SubsysReco
0018
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHIODataNode.h>
0021 #include <phool/PHNode.h>  // for PHNode
0022 #include <phool/PHNodeIterator.h>
0023 #include <phool/PHObject.h>  // for PHObject
0024 #include <phool/getClass.h>
0025 #include <phool/phool.h>  // for PHWHERE
0026
0027 #include <TRotation.h>
0028 #include <TVector3.h>
0029
0030 #include <cstdlib>    // for exit
0031 #include <exception>  // for exception
0032 #include <fstream>
0033 #include <iostream>
0034 #include <map>
0035 #include <sstream>
0036 #include <stdexcept>
0037 #include <utility>  // for pair, make_pair
0038
0039 RawTowerBuilderByHitIndex::RawTowerBuilderByHitIndex(const std::string &name)
0040   : SubsysReco(name)
0041 {
0042 }
0043
0044 int RawTowerBuilderByHitIndex::InitRun(PHCompositeNode *topNode)
0045 {
0046   PHNodeIterator iter(topNode);
0047
0048   // Looking for the DST node
0049   PHCompositeNode *dstNode;
0050   dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0051   if (!dstNode)
0052   {
0053     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0054     exit(1);
0055   }
0056
0057   try
0058   {
0059     CreateNodes(topNode);
0060   }
0061   catch (std::exception &e)
0062   {
0063     std::cout << e.what() << std::endl;
0064     // exit(1);
0065   }
0066
0067   try
0068   {
0069     ReadGeometryFromTable();
0070   }
0071   catch (std::exception &e)
0072   {
0073     std::cout << e.what() << std::endl;
0074     // exit(1);
0075   }
0076
0077   return Fun4AllReturnCodes::EVENT_OK;
0078 }
0079
0080 int RawTowerBuilderByHitIndex::process_event(PHCompositeNode *topNode)
0081 {
0082   // get hits
0083   std::string NodeNameHits = "G4HIT_" + m_Detector;
0084   PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, NodeNameHits);
0085   if (!g4hit)
0086   {
0087     std::cout << "Could not locate g4 hit node " << NodeNameHits << std::endl;
0088     exit(1);
0089   }
0090
0091   // loop over all hits in the event
0092   PHG4HitContainer::ConstIterator hiter;
0093   PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
0094
0095   for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
0096   {
0097     PHG4Hit *g4hit_i = hiter->second;
0098
0099     // Don't include hits with zero energy
0100     if (g4hit_i->get_edep() <= 0 && g4hit_i->get_edep() != -1)
0101     {
0102       continue;
0103     }
0104
0105     /* encode CaloTowerID from j, k index of tower / hit and calorimeter ID */
0106     RawTowerDefs::keytype calotowerid = RawTowerDefs::encode_towerid(m_CaloId,
0107                                                                      g4hit_i->get_index_j(),
0108                                                                      g4hit_i->get_index_k());
0109
0110     /* add the energy to the corresponding tower */
0111     RawTowerv1 *tower = dynamic_cast<RawTowerv1 *>(m_Towers->getTower(calotowerid));
0112     if (!tower)
0113     {
0114       tower = new RawTowerv1(calotowerid);
0115       tower->set_energy(0);
0116       m_Towers->AddTower(tower->get_id(), tower);
0117     }
0118
0119     tower->add_ecell(((unsigned int) (g4hit_i->get_index_j()) << 16U) + g4hit_i->get_index_k(), g4hit_i->get_light_yield());
0120     tower->set_energy(tower->get_energy() + g4hit_i->get_light_yield());
0121     tower->add_eshower(g4hit_i->get_shower_id(), g4hit_i->get_edep());
0122   }
0123
0124   float towerE = 0.;
0125
0126   if (Verbosity())
0127   {
0128     towerE = m_Towers->getTotalEdep();
0129     std::cout << "towers before compression: " << m_Towers->size() << "\t" << m_Detector << std::endl;
0130   }
0131   m_Towers->compress(m_Emin);
0132   if (Verbosity())
0133   {
0134     std::cout << "storing towers: " << m_Towers->size() << std::endl;
0135     std::cout << "Energy lost by dropping towers with less than " << m_Emin
0136               << " energy, lost energy: " << towerE - m_Towers->getTotalEdep() << std::endl;
0137     m_Towers->identify();
0138     RawTowerContainer::ConstRange begin_end = m_Towers->getTowers();
0139     RawTowerContainer::ConstIterator iter;
0140     for (iter = begin_end.first; iter != begin_end.second; ++iter)
0141     {
0142       iter->second->identify();
0143     }
0144   }
0145
0146   return Fun4AllReturnCodes::EVENT_OK;
0147 }
0148
0149 int RawTowerBuilderByHitIndex::End(PHCompositeNode * /*topNode*/)
0150 {
0151   return Fun4AllReturnCodes::EVENT_OK;
0152 }
0153
0154 void RawTowerBuilderByHitIndex::Detector(const std::string &d)
0155 {
0156   m_Detector = d;
0157   m_CaloId = RawTowerDefs::convert_name_to_caloid(m_Detector);
0158 }
0159
0160 void RawTowerBuilderByHitIndex::CreateNodes(PHCompositeNode *topNode)
0161 {
0162   PHNodeIterator iter(topNode);
0163   PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0164   if (!runNode)
0165   {
0166     std::cout << PHWHERE << "Run Node missing, doing nothing." << std::endl;
0167     throw std::runtime_error("Failed to find Run node in RawTowerBuilderByHitIndex::CreateNodes");
0168   }
0169
0170   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0171   if (!dstNode)
0172   {
0173     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0174     throw std::runtime_error("Failed to find DST node in RawTowerBuilderByHitIndex::CreateNodes");
0175   }
0176
0177   // Create the tower geometry node on the tree
0178   m_Geoms = new RawTowerGeomContainerv1(RawTowerDefs::convert_name_to_caloid(m_Detector));
0179   std::string NodeNameTowerGeometries = "TOWERGEOM_" + m_Detector;
0180
0181   PHIODataNode<PHObject> *geomNode = new PHIODataNode<PHObject>(m_Geoms, NodeNameTowerGeometries, "PHObject");
0182   runNode->addNode(geomNode);
0183
0184   // Find detector node (or create new one if not found)
0185   PHNodeIterator dstiter(dstNode);
0186   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst(
0187       "PHCompositeNode", m_Detector));
0188   if (!DetNode)
0189   {
0190     DetNode = new PHCompositeNode(m_Detector);
0191     dstNode->addNode(DetNode);
0192   }
0193
0194   // Create the tower nodes on the tree
0195   m_Towers = new RawTowerContainer(RawTowerDefs::convert_name_to_caloid(m_Detector));
0196   std::string NodeNameTowers;
0197   if (m_SimTowerNodePrefix.empty())
0198   {
0199     // no prefix, consistent with older convension
0200     NodeNameTowers = "TOWER_" + m_Detector;
0201   }
0202   else
0203   {
0204     NodeNameTowers = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
0205   }
0206
0207   PHIODataNode<PHObject> *towerNode = new PHIODataNode<PHObject>(m_Towers, NodeNameTowers, "PHObject");
0208   DetNode->addNode(towerNode);
0209
0210   return;
0211 }
0212
0213 bool RawTowerBuilderByHitIndex::ReadGeometryFromTable()
0214 {
0215   /* Stream to read table from file */
0216   std::ifstream istream_mapping;
0217
0218   /* Open the datafile, if it won't open return an error */
0219   if (!istream_mapping.is_open())
0220   {
0221     istream_mapping.open(m_MappingTowerFile);
0222     if (!istream_mapping)
0223     {
0224       std::cout << "CaloTowerGeomManager::ReadGeometryFromTable - ERROR Failed to open mapping file " << m_MappingTowerFile << std::endl;
0225       exit(1);
0226     }
0227   }
0228
0229   std::string line_mapping;
0230
0231   while (getline(istream_mapping, line_mapping))
0232   {
0233     /* Skip lines starting with / including a '#' */
0234     if (line_mapping.find('#') != std::string::npos)
0235     {
0236       if (Verbosity() > 0)
0237       {
0238         std::cout << "RawTowerBuilderByHitIndex: SKIPPING line in mapping file: " << line_mapping << std::endl;
0239       }
0240       continue;
0241     }
0242
0243     std::istringstream iss(line_mapping);
0244
0245     /* If line starts with keyword Tower, add to tower positions */
0246     if (line_mapping.find("Tower ") != std::string::npos)
0247     {
0248       unsigned idx_j;
0249       unsigned idx_k;
0250       unsigned idx_l;
0251       double pos_x;
0252       double pos_y;
0253       double pos_z;
0254       double size_x;
0255       double size_y;
0256       double size_z;
0257       double rot_x;
0258       double rot_y;
0259       double rot_z;
0260       double type;
0261       std::string dummys;
0262
0263       /* read string- break if error */
0264       if (!(iss >> dummys >> type >> idx_j >> idx_k >> idx_l >> pos_x >> pos_y >> pos_z >> size_x >> size_y >> size_z >> rot_x >> rot_y >> rot_z))
0265       {
0266         std::cout << "ERROR in RawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << std::endl;
0267         exit(1);
0268       }
0269
0270       /* Construct unique Tower ID */
0271       unsigned int temp_id = RawTowerDefs::encode_towerid(m_CaloId, idx_j, idx_k);
0272
0273       /* Create tower geometry object */
0274       RawTowerGeom *temp_geo = new RawTowerGeomv3(temp_id);
0275       temp_geo->set_center_x(pos_x);
0276       temp_geo->set_center_y(pos_y);
0277       temp_geo->set_center_z(pos_z);
0278       temp_geo->set_size_x(size_x);
0279       temp_geo->set_size_y(size_y);
0280       temp_geo->set_size_z(size_z);
0281       temp_geo->set_tower_type((int) type);
0282
0283       /* Insert this tower into position map */
0284       m_Geoms->add_tower_geometry(temp_geo);
0285     }
0286     /* If line does not start with keyword Tower, read as parameter */
0287     else
0288     {
0289       /* If this line is not a comment and not a tower, save parameter as string / value. */
0290       std::string parname;
0291       double parval;
0292
0293       /* read string- break if error */
0294       if (!(iss >> parname >> parval))
0295       {
0296         std::cout << "ERROR in RawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << std::endl;
0297         exit(1);
0298       }
0299
0300       m_GlobalParameterMap.insert(std::make_pair(parname, parval));
0301     }
0302   }
0303
0304   /* Update member variables for global parameters based on parsed parameter file */
0305   std::map<std::string, double>::iterator parit;
0306
0307   parit = m_GlobalParameterMap.find("Gx0");
0308   if (parit != m_GlobalParameterMap.end())
0309   {
0310     m_GlobalPlaceInX = parit->second;
0311   }
0312
0313   parit = m_GlobalParameterMap.find("Gy0");
0314   if (parit != m_GlobalParameterMap.end())
0315   {
0316     m_GlobalPlaceInY = parit->second;
0317   }
0318
0319   parit = m_GlobalParameterMap.find("Gz0");
0320   if (parit != m_GlobalParameterMap.end())
0321   {
0322     m_GlobalPlaceInZ = parit->second;
0323   }
0324
0325   parit = m_GlobalParameterMap.find("Grot_x");
0326   if (parit != m_GlobalParameterMap.end())
0327   {
0328     m_RotInX = parit->second;
0329   }
0330
0331   parit = m_GlobalParameterMap.find("Grot_y");
0332   if (parit != m_GlobalParameterMap.end())
0333   {
0334     m_RotInY = parit->second;
0335   }
0336
0337   parit = m_GlobalParameterMap.find("Grot_z");
0338   if (parit != m_GlobalParameterMap.end())
0339   {
0340     m_RotInZ = parit->second;
0341   }
0342
0343   /* Correct tower geometries for global calorimter translation / rotation
0344    * after reading parameters from file */
0345   RawTowerGeomContainer::ConstRange all_towers = m_Geoms->get_tower_geometries();
0346
0347   for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
0348        it != all_towers.second; ++it)
0349   {
0350     double x_temp = it->second->get_center_x();
0351     double y_temp = it->second->get_center_y();
0352     double z_temp = it->second->get_center_z();
0353
0354     /* Rotation */
0355     TRotation rot;
0356     rot.RotateX(m_RotInX);
0357     rot.RotateY(m_RotInY);
0358     rot.RotateZ(m_RotInZ);
0359
0360     TVector3 v_temp_r(x_temp, y_temp, z_temp);
0361     v_temp_r.Transform(rot);
0362
0363     /* Translation */
0364     double x_temp_rt = v_temp_r.X() + m_GlobalPlaceInX;
0365     double y_temp_rt = v_temp_r.Y() + m_GlobalPlaceInY;
0366     double z_temp_rt = v_temp_r.Z() + m_GlobalPlaceInZ;
0367
0368     /* Update tower geometry object */
0369     it->second->set_center_x(x_temp_rt);
0370     it->second->set_center_y(y_temp_rt);
0371     it->second->set_center_z(z_temp_rt);
0372
0373     if (Verbosity() > 2)
0374     {
0375       std::cout << "* Local tower x y z : " << x_temp << " " << y_temp << " " << z_temp << std::endl;
0376       std::cout << "* Globl tower x y z : " << x_temp_rt << " " << y_temp_rt << " " << z_temp_rt << std::endl;
0377     }
0378   }
0379
0380   if (Verbosity())
0381   {
0382     std::cout << "size tower geom container:" << m_Geoms->size() << "\t" << m_Detector << std::endl;
0383   }
0384   return true;
0385 }