packages/CaloReco/CaloGeomMapping.cc

0001 #include "CaloGeomMapping.h"
0002
0003 #include <cdbobjects/CDBTTree.h>
0004
0005 #include <ffamodules/CDBInterface.h>
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <fun4all/SubsysReco.h>  // for SubsysReco
0009
0010 #include <phool/PHCompositeNode.h>
0011 #include <phool/PHIODataNode.h>  // for PHIODataNode
0012 #include <phool/PHNode.h>
0013 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0014 #include <phool/PHObject.h>
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>  // for PHWHERE
0017
0018 #include <calobase/RawTowerDefs.h>           // for encode_towerid
0019 #include <calobase/RawTowerGeom.h>           // for RawTowerGeom
0020 #include <calobase/RawTowerGeomContainer.h>  // for RawTowerGeomC...
0021 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
0022 #include <calobase/RawTowerGeomv1.h>
0023 #include <calobase/RawTowerGeomv5.h>
0024
0025 #include <cmath>      // for fabs, atan, cos
0026 #include <cstdlib>    // for exit
0027 #include <exception>  // for exception
0028 #include <iostream>   // for operator<<, endl
0029 #include <stdexcept>  // for runtime_error
0030 #include <utility>    // for pair
0031
0032 //____________________________________________________________________________..
0033 CaloGeomMapping::CaloGeomMapping(const std::string &name)
0034   : SubsysReco(name)
0035   , m_Detector("CEMC")
0036 {
0037 }
0038
0039 //____________________________________________________________________________..
0040 int CaloGeomMapping::Init(PHCompositeNode *topNode)
0041 {
0042   if (Verbosity() > 0)
0043   {
0044     std::cout << "CaloGeomMapping::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0045   }
0046   /* std::cout << "Printing node tree before new node creation:" << std::endl; */
0047   /* topNode->print(); */
0048   try
0049   {
0050     CreateGeomNode(topNode);
0051   }
0052   catch (std::exception &e)
0053   {
0054     std::cout << e.what() << std::endl;
0055     exit(1);
0056   }
0057   /* std::cout << "Printing node tree after new node creation:" << std::endl; */
0058   /* topNode->print(); */
0059   return Fun4AllReturnCodes::EVENT_OK;
0060 }
0061
0062 void CaloGeomMapping::CreateGeomNode(PHCompositeNode *topNode)
0063 {
0064   PHNodeIterator iter(topNode);
0065   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0066   if (!runNode)
0067   {
0068     std::cout << PHWHERE << "Run Node missing, doing nothing." << std::endl;
0069     throw std::runtime_error("Failed to find Run node in CaloGeomMapping::CreateGeomNode");
0070   }
0071
0072   PHNodeIterator runIter(runNode);
0073   PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_Detector));
0074   if (!RunDetNode)
0075   {
0076     RunDetNode = new PHCompositeNode(m_Detector);
0077     runNode->addNode(RunDetNode);
0078   }
0079
0080   m_caloid = RawTowerDefs::convert_name_to_caloid(m_Detector);
0081   m_TowerGeomNodeName = "TOWERGEOM_" + m_Detector;
0082
0083   if (m_UseDetailedGeometry && m_Detector == "CEMC")
0084   {
0085     // Rename the TOWERGEOM node to avoid confusion with the former geometry
0086     m_TowerGeomNodeName = m_TowerGeomNodeName + "_DETAILED";
0087   }
0088
0089   m_RawTowerGeomContainer = findNode::getClass<RawTowerGeomContainer>(topNode, m_TowerGeomNodeName);
0090   if (!m_RawTowerGeomContainer)
0091   {
0092     m_RawTowerGeomContainer = new RawTowerGeomContainer_Cylinderv1(m_caloid);
0093     // add it to the node tree
0094     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(m_RawTowerGeomContainer, m_TowerGeomNodeName, "PHObject");
0095     RunDetNode->addNode(newNode);
0096   }
0097
0098   // Get the geometry mapping file from the Conditions Database
0099   std::string inName=CDBInterface::instance()->getUrl("CALO_TOWER_GEOMETRY");
0100   CDBTTree * cdbttree = new CDBTTree(inName);
0101   cdbttree->LoadCalibrations();
0102
0103   std::string parName;
0104   std::string parBase;
0105   // Set the radius, thickness, number of eta and phi bins
0106   if (m_Detector == "CEMC")
0107   {
0108     parBase = "cemc";
0109     m_RawTowerGeomContainer->set_radius(93.5);
0110     m_RawTowerGeomContainer->set_thickness(20.4997);
0111     m_RawTowerGeomContainer->set_phibins(256);
0112     m_RawTowerGeomContainer->set_etabins(96);
0113   }
0114   if (m_Detector == "HCALIN")
0115   {
0116     parBase = "hcalin";
0117     m_RawTowerGeomContainer->set_radius(115);
0118     m_RawTowerGeomContainer->set_thickness(25.005);
0119     m_RawTowerGeomContainer->set_phibins(64);
0120     m_RawTowerGeomContainer->set_etabins(24);
0121   }
0122   if (m_Detector == "HCALOUT")
0123   {
0124     parBase = "hcalout";
0125     m_RawTowerGeomContainer->set_radius(177.423);
0126     m_RawTowerGeomContainer->set_thickness(96.894);
0127     m_RawTowerGeomContainer->set_phibins(64);
0128     m_RawTowerGeomContainer->set_etabins(24);
0129   }
0130
0131   // Set the eta and phi bounds of each bin
0132   for (int ibin = 0; ibin < m_RawTowerGeomContainer->get_etabins(); ibin++)
0133   {
0134     parName = parBase + "_eta_";
0135     double first, second;
0136     first = cdbttree->GetDoubleValue(ibin, parName + "first");
0137     second = cdbttree->GetDoubleValue(ibin, parName + "second");
0138     const std::pair<double, double> range(first, second);
0139     m_RawTowerGeomContainer->set_etabounds(ibin, range);
0140   }
0141   for (int ibin = 0; ibin < m_RawTowerGeomContainer->get_phibins(); ibin++)
0142   {
0143     parName = parBase + "_phi_";
0144     double first, second;
0145     first = cdbttree->GetDoubleValue(ibin, parName + "first");
0146     second = cdbttree->GetDoubleValue(ibin, parName + "second");
0147     const std::pair<double, double> range(first, second);
0148     m_RawTowerGeomContainer->set_phibounds(ibin, range);
0149   }
0150
0151
0152   // Build the RawTowerGeom objects
0153   if (m_UseDetailedGeometry == true && m_Detector == "CEMC")
0154   {
0155     // Detailed geometry using the 8 block vertices as defined in the GEANT4 simulation
0156     BuildDetailedGeometry();
0157   }
0158   else
0159   {
0160     // Approximate geometry
0161     // Deduces the tower position from the eta/phi bins
0162     // and project them at a constant radius.
0163     if (m_UseDetailedGeometry == true)
0164     {
0165       std::cout << "CaloGeomMapping::CreateNodes - Detailed geometry is not yet supported for " << m_Detector << ". The former approximate geometry is used instead." << std::endl;
0166     }
0167     BuildFormerGeometry();
0168   }
0169 }
0170
0171 void CaloGeomMapping::BuildFormerGeometry()
0172 {
0173   // Populate container with RawTowerGeom objects
0174   for (int ieta = 0; ieta < m_RawTowerGeomContainer->get_etabins(); ieta++)
0175   {
0176     for (int iphi = 0; iphi < m_RawTowerGeomContainer->get_phibins(); iphi++)
0177     {
0178       // build tower geom here
0179       const RawTowerDefs::keytype key =
0180           RawTowerDefs::encode_towerid(m_caloid, ieta, iphi);
0181
0182       double r = m_RawTowerGeomContainer->get_radius();
0183       const double x(r * cos(m_RawTowerGeomContainer->get_phicenter(iphi)));
0184       const double y(r * sin(m_RawTowerGeomContainer->get_phicenter(iphi)));
0185       const double z(r / tan(2 * atan(exp(-1 * m_RawTowerGeomContainer->get_etacenter(ieta)))));
0186       RawTowerGeom *tg = m_RawTowerGeomContainer->get_tower_geometry(key);
0187       if (tg)
0188       {
0189         if (Verbosity() > 0)
0190         {
0191           std::cout << "CaloGeomMapping::BuildFormerGeometry - Tower geometry " << key << " already exists" << std::endl;
0192         }
0193
0194         if (fabs(tg->get_center_x() - x) > 1e-4)
0195         {
0196           std::cout << "CaloGeomMapping::BuildFormerGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
0197                     << std::endl;
0198
0199           exit(1);
0200         }
0201         if (fabs(tg->get_center_y() - y) > 1e-4)
0202         {
0203           std::cout << "CaloGeomMapping::BuildFormerGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
0204                     << std::endl;
0205           exit(1);
0206         }
0207         if (fabs(tg->get_center_z() - z) > 1e-4)
0208         {
0209           std::cout << "CaloGeomMapping::BuildFormerGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
0210                     << std::endl;
0211           exit(1);
0212         }
0213       }
0214       else
0215       {
0216         if (Verbosity() > 0)
0217         {
0218           std::cout << "CaloGeomMapping::BuildFormerGeometry - building tower geometry " << key << "" << std::endl;
0219         }
0220
0221         tg = new RawTowerGeomv1(key);
0222
0223         tg->set_center_x(x);
0224         tg->set_center_y(y);
0225         tg->set_center_z(z);
0226         m_RawTowerGeomContainer->add_tower_geometry(tg);
0227       }
0228     }
0229   }  // end loop over eta, phi bins
0230 }
0231
0232 void CaloGeomMapping::BuildDetailedGeometry()
0233 {
0234   // This method is only implemented for the electromagnetic calorimeter
0235
0236   // Get the geometry mapping file from the Conditions Database
0237   std::string inName = CDBInterface::instance()->getUrl("CALO_TOWER_GEOMETRY_EMCAL_DETAILED");
0238   CDBTTree *cdbttree = new CDBTTree(inName);
0239   cdbttree->LoadCalibrations();
0240   //std::cout << "printing CDB TTree" << std::endl;
0241   //cdbttree->Print();
0242
0243   // Populate container with RawTowerGeom objects
0244   const int nSectors = 32;
0245   const int nUniqueTowers = 192;
0246   for (int iSector = 0; iSector < nSectors; iSector++)
0247   {
0248     for (int iUniqueTower = 0; iUniqueTower < nUniqueTowers; iUniqueTower++)
0249     {
0250       // Vertices' coordinates of the corresponding unique block
0251       // separated into 4 towers (2 bins in eta x 2 bins in phi)
0252
0253       /*   R
0254            |                    7_______14      14_______6
0255            |                   /|      /|       /|      /|
0256            |_______ eta     15/_|_____/18    18/_|_____/13
0257           /                   | |     | |      | |     | |
0258          /                    |3|_____|_|10    10|_____|_|2
0259         /                     | /     | /      | /     | /
0260        phi                  11|/______|/17   17|/______|/9
0261
0262                           15_______18      18_______13
0263                           /|      /|       /|      /|
0264                         8/_|_____/16    16/_|_____/5|
0265                          | |     | |      | |     | |
0266                          11|_____|_|17    17|_____|_|9
0267                          | /     | /      | /     | /
0268                         4|/______|/12   12|/______|/1
0269       */
0270
0271
0272       std::vector<double> vertices(18*3);
0273       // Block vertices (1 -> 8)
0274       for (int i = 0; i < 8; i++)
0275       {
0276         vertices[i*3+0] = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1, "vtx_" + std::to_string(i) + "_x") / 10; // Convert mm to cm
0277         vertices[i*3+1] = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1, "vtx_" + std::to_string(i) + "_y") / 10;
0278         vertices[i*3+2] = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1, "vtx_" + std::to_string(i) + "_z") / 10;
0279       }
0280
0281       // Get the rotation of the tower:
0282       double rot_x = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1,"rot_x");
0283       double rot_y = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1,"rot_y");
0284       double rot_z = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1,"rot_z");
0285
0286       // Divide each block into 4 equal towers.
0287       for (int j = 0; j < 3; j++)
0288       {
0289         for (int i = 0; i < 4; i++)
0290         {
0291           vertices[(8 + i) * 3 + j]=(vertices[i * 3 + j] + vertices[((i + 1) % 4) * 3 + j]) / 2; // Middle points corresponding to tower vertices (inner side)
0292           vertices[(8 + i + 4) * 3 + j]=(vertices[(i + 4) * 3 + j] + vertices[((i + 1) % 4 + 4) * 3 + j]) / 2; // Middle points (outer side)
0293         }
0294         vertices[(16)*3+j] = (vertices[0 * 3 + j] + vertices[1 * 3 + j] + vertices[2 * 3 + j] + vertices[3 * 3 + j]) / 4;
0295         vertices[(17)*3+j] = (vertices[4 * 3 + j] + vertices[5 * 3 + j] + vertices[6 * 3 + j] + vertices[7 * 3 + j]) / 4;
0296       }
0297
0298       int ietaBlock = iUniqueTower / 4;
0299       int iphiBlock = iSector * 4 + iUniqueTower % 4;
0300
0301       // Vertices' indices for each of the 4 individual towers.
0302       double sub_tower[2][2][8];
0303       double sub_tower_01[8] = {12, 17, 11, 4, 16, 18, 15, 8};
0304       double sub_tower_00[8] = {17, 10, 3, 11, 18, 14, 7, 15};
0305       double sub_tower_11[8] = {1, 9, 17, 12, 5, 13, 18, 16};
0306       double sub_tower_10[8] = {9, 2, 10, 17, 13, 6, 14, 18};
0307       for (int ivtx = 0; ivtx < 8; ivtx++)
0308       {
0309         sub_tower[0][0][ivtx] = sub_tower_00[ivtx];
0310         sub_tower[0][1][ivtx] = sub_tower_01[ivtx];
0311         sub_tower[1][0][ivtx] = sub_tower_10[ivtx];
0312         sub_tower[1][1][ivtx] = sub_tower_11[ivtx];
0313       }
0314
0315       for (int iTower = 0; iTower < 2; iTower++)
0316       {
0317         for (int jTower = 0; jTower < 2; jTower++)
0318         {
0319           int ieta = ietaBlock * 2 + iTower;
0320           int iphi = iphiBlock * 2 + jTower;
0321
0322           std::vector<double> towerVertices(8*3);
0323           for (int ivtx = 0; ivtx < 8; ivtx++)
0324           {
0325             for (int icoord = 0; icoord < 3; icoord++)
0326             {
0327               towerVertices[ivtx*3+icoord] = vertices[(sub_tower[iTower][jTower][ivtx]-1)*3+icoord];
0328             }
0329           }
0330
0331           // build tower geom here
0332           const RawTowerDefs::keytype key =
0333             RawTowerDefs::encode_towerid(m_caloid, ieta, iphi);
0334
0335
0336           RawTowerGeomv5 *tg0 = new RawTowerGeomv5(key);
0337
0338           tg0->set_vertices(towerVertices);
0339           tg0->set_rotx(rot_x);
0340           tg0->set_roty(rot_y);
0341           tg0->set_rotz(rot_z);
0342
0343           const double x(tg0->get_center_x());
0344           const double y(tg0->get_center_y());
0345           const double z(tg0->get_center_z());
0346
0347           RawTowerGeom *tg = m_RawTowerGeomContainer->get_tower_geometry(key);
0348           if (tg)
0349           {
0350             delete tg0;
0351             if (Verbosity() > 0)
0352             {
0353               std::cout << "CaloGeomMapping::BuildDetailedGeometry - Tower geometry " << key << " already exists" << std::endl;
0354             }
0355
0356             if (fabs(tg->get_center_x() - x) > 1e-4)
0357             {
0358               std::cout << "CaloGeomMapping::BuildDetailedGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
0359                         << std::endl;
0360
0361               exit(1);
0362             }
0363             if (fabs(tg->get_center_y() - y) > 1e-4)
0364             {
0365               std::cout << "CaloGeomMapping::BuildDetailedGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
0366                         << std::endl;
0367               exit(1);
0368             }
0369             if (fabs(tg->get_center_z() - z) > 1e-4)
0370             {
0371               std::cout << "CaloGeomMapping::BuildDetailedGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
0372                         << std::endl;
0373               exit(1);
0374             }
0375           }
0376           else
0377           {
0378             if (Verbosity() > 0)
0379             {
0380               std::cout << "CaloGeomMapping::BuildDetailedGeometry - building tower geometry " << key << "" << std::endl;
0381             }
0382
0383             tg = tg0;
0384             m_RawTowerGeomContainer->add_tower_geometry(tg);
0385           }
0386         }
0387       }
0388     }
0389   }  // end loop over eta, phi bins
0390 }
0391