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;
0136     double second;
0137     first = cdbttree->GetDoubleValue(ibin, parName + "first");
0138     second = cdbttree->GetDoubleValue(ibin, parName + "second");
0139     const std::pair<double, double> range(first, second);
0140     m_RawTowerGeomContainer->set_etabounds(ibin, range);
0141   }
0142   for (int ibin = 0; ibin < m_RawTowerGeomContainer->get_phibins(); ibin++)
0143   {
0144     parName = parBase + "_phi_";
0145     double first;
0146     double second;
0147     first = cdbttree->GetDoubleValue(ibin, parName + "first");
0148     second = cdbttree->GetDoubleValue(ibin, parName + "second");
0149     const std::pair<double, double> range(first, second);
0150     m_RawTowerGeomContainer->set_phibounds(ibin, range);
0151   }
0152
0153
0154   // Build the RawTowerGeom objects
0155   if (m_UseDetailedGeometry == true && m_Detector == "CEMC")
0156   {
0157     // Detailed geometry using the 8 block vertices as defined in the GEANT4 simulation
0158     BuildDetailedGeometry();
0159   }
0160   else
0161   {
0162     // Approximate geometry
0163     // Deduces the tower position from the eta/phi bins
0164     // and project them at a constant radius.
0165     if (m_UseDetailedGeometry == true)
0166     {
0167       std::cout << "CaloGeomMapping::CreateNodes - Detailed geometry is not yet supported for " << m_Detector << ". The former approximate geometry is used instead." << std::endl;
0168     }
0169     BuildFormerGeometry();
0170   }
0171 }
0172
0173 void CaloGeomMapping::BuildFormerGeometry()
0174 {
0175   // Populate container with RawTowerGeom objects
0176   for (int ieta = 0; ieta < m_RawTowerGeomContainer->get_etabins(); ieta++)
0177   {
0178     for (int iphi = 0; iphi < m_RawTowerGeomContainer->get_phibins(); iphi++)
0179     {
0180       // build tower geom here
0181       const RawTowerDefs::keytype key =
0182           RawTowerDefs::encode_towerid(m_caloid, ieta, iphi);
0183
0184       double r = m_RawTowerGeomContainer->get_radius();
0185
0186       if (m_Detector == "HCALIN" || m_Detector == "HCALOUT")
0187       {
0188         r += m_RawTowerGeomContainer->get_radius() / 2.;
0189       }
0190
0191       const double x(r * cos(m_RawTowerGeomContainer->get_phicenter(iphi)));
0192       const double y(r * sin(m_RawTowerGeomContainer->get_phicenter(iphi)));
0193       const double z(r / tan(2 * atan(exp(-1 * m_RawTowerGeomContainer->get_etacenter(ieta)))));
0194       RawTowerGeom *tg = m_RawTowerGeomContainer->get_tower_geometry(key);
0195       if (tg)
0196       {
0197         if (Verbosity() > 0)
0198         {
0199           std::cout << "CaloGeomMapping::BuildFormerGeometry - Tower geometry " << key << " already exists" << std::endl;
0200         }
0201
0202         if (fabs(tg->get_center_x() - x) > 1e-4)
0203         {
0204           std::cout << "CaloGeomMapping::BuildFormerGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
0205                     << std::endl;
0206
0207           exit(1);
0208         }
0209         if (fabs(tg->get_center_y() - y) > 1e-4)
0210         {
0211           std::cout << "CaloGeomMapping::BuildFormerGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
0212                     << std::endl;
0213           exit(1);
0214         }
0215         if (fabs(tg->get_center_z() - z) > 1e-4)
0216         {
0217           std::cout << "CaloGeomMapping::BuildFormerGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
0218                     << std::endl;
0219           exit(1);
0220         }
0221       }
0222       else
0223       {
0224         if (Verbosity() > 0)
0225         {
0226           std::cout << "CaloGeomMapping::BuildFormerGeometry - building tower geometry " << key << "" << std::endl;
0227         }
0228
0229         tg = new RawTowerGeomv1(key);
0230
0231         tg->set_center_x(x);
0232         tg->set_center_y(y);
0233         tg->set_center_z(z);
0234         m_RawTowerGeomContainer->add_tower_geometry(tg);
0235       }
0236     }
0237   }  // end loop over eta, phi bins
0238 }
0239
0240 void CaloGeomMapping::BuildDetailedGeometry()
0241 {
0242   // This method is only implemented for the electromagnetic calorimeter
0243
0244   // Get the geometry mapping file from the Conditions Database
0245   std::string inName = CDBInterface::instance()->getUrl("CALO_TOWER_GEOMETRY_EMCAL_DETAILED");
0246   CDBTTree *cdbttree = new CDBTTree(inName);
0247   cdbttree->LoadCalibrations();
0248   //std::cout << "printing CDB TTree" << std::endl;
0249   //cdbttree->Print();
0250
0251   // Populate container with RawTowerGeom objects
0252   const int nSectors = 32;
0253   const int nUniqueTowers = 192;
0254   for (int iSector = 0; iSector < nSectors; iSector++)
0255   {
0256     for (int iUniqueTower = 0; iUniqueTower < nUniqueTowers; iUniqueTower++)
0257     {
0258       // Vertices' coordinates of the corresponding unique block
0259       // separated into 4 towers (2 bins in eta x 2 bins in phi)
0260
0261       /*   R
0262            |                    7_______14      14_______6
0263            |                   /|      /|       /|      /|
0264            |_______ eta     15/_|_____/18    18/_|_____/13
0265           /                   | |     | |      | |     | |
0266          /                    |3|_____|_|10    10|_____|_|2
0267         /                     | /     | /      | /     | /
0268        phi                  11|/______|/17   17|/______|/9
0269
0270                           15_______18      18_______13
0271                           /|      /|       /|      /|
0272                         8/_|_____/16    16/_|_____/5|
0273                          | |     | |      | |     | |
0274                          11|_____|_|17    17|_____|_|9
0275                          | /     | /      | /     | /
0276                         4|/______|/12   12|/______|/1
0277       */
0278
0279
0280       std::vector<double> vertices(18*3);
0281       // Block vertices (1 -> 8)
0282       for (int i = 0; i < 8; i++)
0283       {
0284         vertices[i*3+0] = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1, "vtx_" + std::to_string(i) + "_x") / 10; // Convert mm to cm
0285         vertices[i*3+1] = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1, "vtx_" + std::to_string(i) + "_y") / 10;
0286         vertices[i*3+2] = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1, "vtx_" + std::to_string(i) + "_z") / 10;
0287       }
0288
0289       // Get the rotation of the tower:
0290       double rot_x = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1,"rot_x");
0291       double rot_y = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1,"rot_y");
0292       double rot_z = cdbttree->GetDoubleValue(iUniqueTower * nSectors + iSector + 1,"rot_z");
0293
0294       // Divide each block into 4 equal towers.
0295       for (int j = 0; j < 3; j++)
0296       {
0297         for (int i = 0; i < 4; i++)
0298         {
0299           vertices[(8 + i) * 3 + j]=(vertices[i * 3 + j] + vertices[((i + 1) % 4) * 3 + j]) / 2; // Middle points corresponding to tower vertices (inner side)
0300           vertices[(8 + i + 4) * 3 + j]=(vertices[(i + 4) * 3 + j] + vertices[((i + 1) % 4 + 4) * 3 + j]) / 2; // Middle points (outer side)
0301         }
0302         vertices[(16)*3+j] = (vertices[0 * 3 + j] + vertices[1 * 3 + j] + vertices[2 * 3 + j] + vertices[3 * 3 + j]) / 4;
0303         vertices[(17)*3+j] = (vertices[4 * 3 + j] + vertices[5 * 3 + j] + vertices[6 * 3 + j] + vertices[7 * 3 + j]) / 4;
0304       }
0305
0306       int ietaBlock = iUniqueTower / 4;
0307       int iphiBlock = iSector * 4 + iUniqueTower % 4;
0308
0309       // Vertices' indices for each of the 4 individual towers.
0310       double sub_tower[2][2][8];
0311       double sub_tower_01[8] = {12, 17, 11, 4, 16, 18, 15, 8};
0312       double sub_tower_00[8] = {17, 10, 3, 11, 18, 14, 7, 15};
0313       double sub_tower_11[8] = {1, 9, 17, 12, 5, 13, 18, 16};
0314       double sub_tower_10[8] = {9, 2, 10, 17, 13, 6, 14, 18};
0315       for (int ivtx = 0; ivtx < 8; ivtx++)
0316       {
0317         sub_tower[0][0][ivtx] = sub_tower_00[ivtx];
0318         sub_tower[0][1][ivtx] = sub_tower_01[ivtx];
0319         sub_tower[1][0][ivtx] = sub_tower_10[ivtx];
0320         sub_tower[1][1][ivtx] = sub_tower_11[ivtx];
0321       }
0322
0323       for (int iTower = 0; iTower < 2; iTower++)
0324       {
0325         for (int jTower = 0; jTower < 2; jTower++)
0326         {
0327           int ieta = ietaBlock * 2 + iTower;
0328           int iphi = iphiBlock * 2 + jTower;
0329
0330           std::vector<double> towerVertices(8*3);
0331           for (int ivtx = 0; ivtx < 8; ivtx++)
0332           {
0333             for (int icoord = 0; icoord < 3; icoord++)
0334             {
0335               towerVertices[ivtx*3+icoord] = vertices[(sub_tower[iTower][jTower][ivtx]-1)*3+icoord];
0336             }
0337           }
0338
0339           // build tower geom here
0340           const RawTowerDefs::keytype key =
0341             RawTowerDefs::encode_towerid(m_caloid, ieta, iphi);
0342
0343
0344           RawTowerGeomv5 *tg0 = new RawTowerGeomv5(key);
0345
0346           tg0->set_vertices(towerVertices);
0347           tg0->set_rotx(rot_x);
0348           tg0->set_roty(rot_y);
0349           tg0->set_rotz(rot_z);
0350
0351           const double x(tg0->get_center_x());
0352           const double y(tg0->get_center_y());
0353           const double z(tg0->get_center_z());
0354
0355           RawTowerGeom *tg = m_RawTowerGeomContainer->get_tower_geometry(key);
0356           if (tg)
0357           {
0358             delete tg0;
0359             if (Verbosity() > 0)
0360             {
0361               std::cout << "CaloGeomMapping::BuildDetailedGeometry - Tower geometry " << key << " already exists" << std::endl;
0362             }
0363
0364             if (fabs(tg->get_center_x() - x) > 1e-4)
0365             {
0366               std::cout << "CaloGeomMapping::BuildDetailedGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
0367                         << std::endl;
0368
0369               exit(1);
0370             }
0371             if (fabs(tg->get_center_y() - y) > 1e-4)
0372             {
0373               std::cout << "CaloGeomMapping::BuildDetailedGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
0374                         << std::endl;
0375               exit(1);
0376             }
0377             if (fabs(tg->get_center_z() - z) > 1e-4)
0378             {
0379               std::cout << "CaloGeomMapping::BuildDetailedGeometry - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
0380                         << std::endl;
0381               exit(1);
0382             }
0383           }
0384           else
0385           {
0386             if (Verbosity() > 0)
0387             {
0388               std::cout << "CaloGeomMapping::BuildDetailedGeometry - building tower geometry " << key << "" << std::endl;
0389             }
0390
0391             tg = tg0;
0392             m_RawTowerGeomContainer->add_tower_geometry(tg);
0393           }
0394         }
0395       }
0396     }
0397   }  // end loop over eta, phi bins
0398 }
0399