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, idx_k, idx_l;
0249       double pos_x, pos_y, pos_z;
0250       double size_x, size_y, size_z;
0251       double rot_x, rot_y, rot_z;
0252       double type;
0253       std::string dummys;
0254
0255       /* read string- break if error */
0256       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))
0257       {
0258         std::cout << "ERROR in RawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << std::endl;
0259         exit(1);
0260       }
0261
0262       /* Construct unique Tower ID */
0263       unsigned int temp_id = RawTowerDefs::encode_towerid(m_CaloId, idx_j, idx_k);
0264
0265       /* Create tower geometry object */
0266       RawTowerGeom *temp_geo = new RawTowerGeomv3(temp_id);
0267       temp_geo->set_center_x(pos_x);
0268       temp_geo->set_center_y(pos_y);
0269       temp_geo->set_center_z(pos_z);
0270       temp_geo->set_size_x(size_x);
0271       temp_geo->set_size_y(size_y);
0272       temp_geo->set_size_z(size_z);
0273       temp_geo->set_tower_type((int) type);
0274
0275       /* Insert this tower into position map */
0276       m_Geoms->add_tower_geometry(temp_geo);
0277     }
0278     /* If line does not start with keyword Tower, read as parameter */
0279     else
0280     {
0281       /* If this line is not a comment and not a tower, save parameter as string / value. */
0282       std::string parname;
0283       double parval;
0284
0285       /* read string- break if error */
0286       if (!(iss >> parname >> parval))
0287       {
0288         std::cout << "ERROR in RawTowerBuilderByHitIndex: Failed to read line in mapping file " << m_MappingTowerFile << std::endl;
0289         exit(1);
0290       }
0291
0292       m_GlobalParameterMap.insert(std::make_pair(parname, parval));
0293     }
0294   }
0295
0296   /* Update member variables for global parameters based on parsed parameter file */
0297   std::map<std::string, double>::iterator parit;
0298
0299   parit = m_GlobalParameterMap.find("Gx0");
0300   if (parit != m_GlobalParameterMap.end())
0301   {
0302     m_GlobalPlaceInX = parit->second;
0303   }
0304
0305   parit = m_GlobalParameterMap.find("Gy0");
0306   if (parit != m_GlobalParameterMap.end())
0307   {
0308     m_GlobalPlaceInY = parit->second;
0309   }
0310
0311   parit = m_GlobalParameterMap.find("Gz0");
0312   if (parit != m_GlobalParameterMap.end())
0313   {
0314     m_GlobalPlaceInZ = parit->second;
0315   }
0316
0317   parit = m_GlobalParameterMap.find("Grot_x");
0318   if (parit != m_GlobalParameterMap.end())
0319   {
0320     m_RotInX = parit->second;
0321   }
0322
0323   parit = m_GlobalParameterMap.find("Grot_y");
0324   if (parit != m_GlobalParameterMap.end())
0325   {
0326     m_RotInY = parit->second;
0327   }
0328
0329   parit = m_GlobalParameterMap.find("Grot_z");
0330   if (parit != m_GlobalParameterMap.end())
0331   {
0332     m_RotInZ = parit->second;
0333   }
0334
0335   /* Correct tower geometries for global calorimter translation / rotation
0336    * after reading parameters from file */
0337   RawTowerGeomContainer::ConstRange all_towers = m_Geoms->get_tower_geometries();
0338
0339   for (RawTowerGeomContainer::ConstIterator it = all_towers.first;
0340        it != all_towers.second; ++it)
0341   {
0342     double x_temp = it->second->get_center_x();
0343     double y_temp = it->second->get_center_y();
0344     double z_temp = it->second->get_center_z();
0345
0346     /* Rotation */
0347     TRotation rot;
0348     rot.RotateX(m_RotInX);
0349     rot.RotateY(m_RotInY);
0350     rot.RotateZ(m_RotInZ);
0351
0352     TVector3 v_temp_r(x_temp, y_temp, z_temp);
0353     v_temp_r.Transform(rot);
0354
0355     /* Translation */
0356     double x_temp_rt = v_temp_r.X() + m_GlobalPlaceInX;
0357     double y_temp_rt = v_temp_r.Y() + m_GlobalPlaceInY;
0358     double z_temp_rt = v_temp_r.Z() + m_GlobalPlaceInZ;
0359
0360     /* Update tower geometry object */
0361     it->second->set_center_x(x_temp_rt);
0362     it->second->set_center_y(y_temp_rt);
0363     it->second->set_center_z(z_temp_rt);
0364
0365     if (Verbosity() > 2)
0366     {
0367       std::cout << "* Local tower x y z : " << x_temp << " " << y_temp << " " << z_temp << std::endl;
0368       std::cout << "* Globl tower x y z : " << x_temp_rt << " " << y_temp_rt << " " << z_temp_rt << std::endl;
0369     }
0370   }
0371
0372   if (Verbosity())
0373   {
0374     std::cout << "size tower geom container:" << m_Geoms->size() << "\t" << m_Detector << std::endl;
0375   }
0376   return true;
0377 }