Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:19:59 

0001 #include "RawClusterBuilderTemplate.h"
0002 
0003 #include "BEmcCluster.h"
0004 #include "BEmcRec.h"
0005 #include "BEmcRecCEMC.h"
0006 
0007 #include <globalvertex/GlobalVertexMap.h>
0008 #include <globalvertex/MbdVertex.h>
0009 #include <globalvertex/MbdVertexMap.h>
0010 // To isolate the issue with the vertex
0011 #include <g4main/PHG4VtxPoint.h>
0012 #include <g4main/PHG4TruthInfoContainer.h>
0013 
0014 #include <calobase/RawCluster.h>
0015 #include <calobase/RawClusterContainer.h>
0016 #include <calobase/RawClusterv1.h>
0017 #include <calobase/RawTower.h>
0018 #include <calobase/RawTowerContainer.h>
0019 #include <calobase/RawTowerDefs.h>
0020 #include <calobase/RawTowerGeom.h>
0021 #include <calobase/RawTowerGeomContainer.h>
0022 #include <calobase/TowerInfo.h>
0023 #include <calobase/TowerInfoContainer.h>
0024 
0025 
0026 #include <ffamodules/CDBInterface.h>
0027 
0028 #include <fun4all/Fun4AllReturnCodes.h>
0029 #include <fun4all/SubsysReco.h>
0030 
0031 #include <phool/PHCompositeNode.h>
0032 #include <phool/PHIODataNode.h>
0033 #include <phool/PHNode.h>
0034 #include <phool/PHNodeIterator.h>
0035 #include <phool/PHObject.h>
0036 #include <phool/getClass.h>
0037 #include <phool/phool.h>
0038 
0039 #include <algorithm>
0040 #include <cmath>
0041 #include <exception>
0042 #include <fstream>
0043 #include <iostream>
0044 #include <map>
0045 #include <stdexcept>
0046 #include <utility>
0047 #include <vector>
0048 #include <limits>
0049 
0050 RawClusterBuilderTemplate::RawClusterBuilderTemplate(const std::string &name)
0051   : SubsysReco(name)
0052 {
0053 }
0054 
0055 RawClusterBuilderTemplate::~RawClusterBuilderTemplate()
0056 {
0057   // one can delete null pointers
0058   delete bemc;
0059 }
0060 
0061 void RawClusterBuilderTemplate::Detector(const std::string &d)
0062 {
0063   detector = d;
0064 
0065   // Create proper BEmcRec object
0066 
0067   if (detector == "CEMC")
0068   {
0069     bemc = new BEmcRecCEMC();
0070   }
0071   else
0072   {
0073     std::cout << "Warning from RawClusterBuilderTemplate::Detector(): no detector specific class "
0074               << Name() << " defined for detector " << detector
0075               << ". Default BEmcRec will be used" << std::endl;
0076     bemc = new BEmcRec();
0077   }
0078 
0079   // Set vertex
0080   float vertex[3] = {0, 0, 0};
0081   bemc->SetVertex(vertex);
0082   // Set threshold
0083   bemc->SetTowerThreshold(_min_tower_e);
0084   bemc->SetPeakThreshold(_min_peak_e);
0085   bemc->SetProbNoiseParam(fProbNoiseParam);
0086 }
0087 
0088 void RawClusterBuilderTemplate::set_UseCorrectPosition(const bool useCorrectPosition)
0089 {
0090   if (bemc == nullptr)
0091   {
0092     std::cerr << "Error in RawClusterBuilderTemplate::set_UseCorrectPosition()(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << std::endl;
0093     return;
0094   }
0095 
0096   bemc->set_UseCorrectPosition(useCorrectPosition);
0097 }
0098 
0099 void RawClusterBuilderTemplate::set_UseCorrectShowerDepth(const bool useCorrectShowerDepth) {
0100   if (bemc == nullptr)
0101   {
0102     std::cerr << "Error in RawClusterBuilderTemplate::set_UseCorrectShowerDepth()(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << std::endl;
0103     return;
0104   }
0105 
0106   bemc->set_UseCorrectShowerDepth(useCorrectShowerDepth);
0107 }
0108 
0109 void RawClusterBuilderTemplate::set_UseDetailedGeometry(const bool useDetailedGeometry)
0110 {
0111   if (bemc == nullptr)
0112   {
0113     std::cerr << "Error in RawClusterBuilderTemplate::set_UseDetailedGeometry()(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << std::endl;
0114     return;
0115   }
0116 
0117   m_UseDetailedGeometry = useDetailedGeometry;
0118   bemc->set_UseDetailedGeometry(m_UseDetailedGeometry);
0119 }
0120 
0121 void RawClusterBuilderTemplate::LoadProfile(const std::string &fname)
0122 {
0123   if (bemc == nullptr)
0124   {
0125     std::cerr << "Error in RawClusterBuilderTemplate::LoadProfile()(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << std::endl;
0126     return;
0127   }
0128   std::string url = CDBInterface::instance()->getUrl("EMCPROFILE", fname);
0129   bemc->LoadProfile(url);
0130 }
0131 
0132 void RawClusterBuilderTemplate::SetCylindricalGeometry()
0133 {
0134   if (bemc == nullptr)
0135   {
0136     std::cerr << "Error in RawClusterBuilderTemplate::SetCylindricalGeometry()(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << std::endl;
0137     return;
0138   }
0139 
0140   bemc->SetCylindricalGeometry();
0141 }
0142 
0143 void RawClusterBuilderTemplate::SetPlanarGeometry()
0144 {
0145   if (bemc == nullptr)
0146   {
0147     std::cerr << "Error in RawClusterBuilderTemplate::SetPlanarGeometry()(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << std::endl;
0148     return;
0149   }
0150 
0151   bemc->SetPlanarGeometry();
0152 }
0153 
0154 int RawClusterBuilderTemplate::InitRun(PHCompositeNode *topNode)
0155 {
0156   if (bemc == nullptr)
0157   {
0158     std::cerr << "Error in RawClusterBuilderTemplate::InitRun(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << std::endl;
0159     return Fun4AllReturnCodes::ABORTEVENT;
0160   }
0161 
0162   // Ensure that the detailed geometry is available if the user requests it.
0163   // Otherwise, use the default geometry 
0164   if (m_UseDetailedGeometry && detector != "CEMC")
0165   {
0166     m_UseDetailedGeometry = false;
0167     bemc->set_UseDetailedGeometry(false);
0168     std::cout << "Warning in RawClusterBuilderTemplate::InitRun()(): No alternative detailed geometry defined for detector " << detector << ". m_UseDetailedGeometry automatically set to false." << std::endl;
0169   }
0170 
0171   try
0172   {
0173     CreateNodes(topNode);
0174   }
0175   catch (std::exception &e)
0176   {
0177     std::cout << PHWHERE << ": " << e.what() << std::endl;
0178     throw;
0179   }
0180 
0181   // geometry case
0182   if (m_TowerGeomNodeName.empty())
0183   {
0184     m_TowerGeomNodeName = "TOWERGEOM_" + detector;
0185     if (m_UseDetailedGeometry)
0186     {
0187       if (detector == "CEMC")
0188       {
0189         m_TowerGeomNodeName = m_TowerGeomNodeName + "_DETAILED";
0190       }
0191       else
0192       {
0193         std::cout << "RawClusterBuilderTemplate::InitRun - Detailed geometry not implemented for detector " << detector << ". The former geometry is used instead" << std::endl;
0194         m_UseDetailedGeometry = false;
0195         bemc->set_UseDetailedGeometry(false);
0196       }
0197     }
0198   }
0199   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, m_TowerGeomNodeName);
0200 
0201   if (!towergeom && m_UseDetailedGeometry)
0202   {
0203     std::cout << "RawClusterBuilderTemplate::InitRun - Detailed geometry node " << m_TowerGeomNodeName << " is not available. "
0204               << "Switching to the former geometry node TOWERGEOM_" << detector << "." << std::endl;
0205     m_UseDetailedGeometry = false;
0206     m_TowerGeomNodeName = "TOWERGEOM_" + detector;
0207     bemc->set_UseDetailedGeometry(false);
0208     towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, m_TowerGeomNodeName);
0209   }
0210   
0211   if (!towergeom)
0212   {
0213     std::cout << PHWHERE << ": Could not find node " << m_TowerGeomNodeName << std::endl;
0214     return Fun4AllReturnCodes::ABORTEVENT;
0215   }
0216 
0217   int Calo_ID = towergeom->get_calorimeter_id();
0218   // std::cout << std::endl << std::endl << std::endl << "Calorimeter ID: " << Calo_ID << std::endl << std::endl << std::endl;
0219 
0220   int ngeom = 0;
0221   int ixmin = 999999;
0222   int ixmax = -999999;
0223   int iymin = 999999;
0224   int iymax = -999999;
0225   RawTowerGeomContainer::ConstRange begin_end_geom = towergeom->get_tower_geometries();
0226   RawTowerGeomContainer::ConstIterator itr_geom = begin_end_geom.first;
0227   for (; itr_geom != begin_end_geom.second; ++itr_geom)
0228   {
0229     RawTowerGeom *towerg = itr_geom->second;
0230     RawTowerDefs::keytype towerid = towerg->get_id();
0231     int ix = RawTowerDefs::decode_index2(towerid);  // index2 is phi in CYL
0232     int iy = RawTowerDefs::decode_index1(towerid);  // index1 is eta in CYL
0233     ixmin = std::min(ixmin, ix);
0234     ixmax = std::max(ixmax, ix);
0235     iymin = std::min(iymin, iy);
0236     iymax = std::max(iymax, iy);
0237     ngeom++;
0238   }
0239   if (Verbosity() > 1)
0240   {
0241     std::cout << "Info from RawClusterBuilderTemplate::InitRun(): Init geometry for "
0242               << detector << ": N of geom towers: " << ngeom << "; ix = "
0243               << ixmin << "-" << ixmax << ", iy = "
0244               << iymin << "-" << iymax << std::endl;
0245   }
0246   if (ixmax < ixmin || iymax < iymin)
0247   {
0248     std::cout << "Error in RawClusterBuilderTemplate::InitRun(): wrong geometry data for detector "
0249               << detector << std::endl;
0250     return Fun4AllReturnCodes::ABORTEVENT;
0251   }
0252 
0253   BINX0 = ixmin;
0254   NBINX = ixmax - ixmin + 1;
0255   BINY0 = iymin;
0256   NBINY = iymax - iymin + 1;
0257 
0258   bemc->SetDim(NBINX, NBINY);
0259 
0260   itr_geom = begin_end_geom.first;
0261   for (; itr_geom != begin_end_geom.second; ++itr_geom)
0262   {
0263     RawTowerGeom *towerg = itr_geom->second;
0264     RawTowerDefs::keytype towerid = towerg->get_id();
0265     //    int itype = towerg->get_tower_type();
0266     //    if( itype==2 ) { // PbSc
0267     int ix = RawTowerDefs::decode_index2(towerid);  // index2 is phi in CYL
0268     int iy = RawTowerDefs::decode_index1(towerid);  // index1 is eta in CYL
0269     ix -= BINX0;
0270     iy -= BINY0;
0271 
0272     // reconstruction change
0273     bemc->SetTowerGeometry(ix, iy, *towerg);
0274     bemc->SetCalotype(Calo_ID);
0275     if (Calo_ID == RawTowerDefs::EEMC ||
0276         Calo_ID == RawTowerDefs::EEMC_crystal ||
0277         Calo_ID == RawTowerDefs::EEMC_glass)
0278     {
0279       bemc->SetScinSize(towerg->get_size_z());
0280     }
0281   }
0282 
0283 
0284   // geometry case
0285   // With the former geometry, the tower dimensions are approximated from
0286   // the consecutive tower center postions in the "CompleteTowerGeometry" method.
0287   // With the detailed geometry, the information is already given in the 
0288   // RawTowerGeom node so no further step is required.
0289   if (!(m_UseDetailedGeometry && detector == "CEMC"))
0290   {
0291     if (!bemc->CompleteTowerGeometry())
0292     {
0293       return Fun4AllReturnCodes::ABORTEVENT;
0294     }
0295   }
0296 
0297   // geometry case
0298   if (bPrintGeom)
0299   {
0300     std::string fname;
0301     if (m_UseDetailedGeometry && detector == "CEMC")
0302     {
0303       fname = "geom_" + detector + "_detailed.txt";
0304       bemc->PrintTowerGeometryDetailed(fname);
0305     }
0306     else
0307     {
0308       fname = "geom_" + detector + ".txt";
0309       bemc->PrintTowerGeometry(fname);
0310     }
0311   }
0312 
0313   // Release memory taken by the RawTowerGeom objects in BEmcRec
0314   // Does nothing if the former geometry is used
0315   bemc->ClearInitialDetailedGeometry();
0316   
0317   return Fun4AllReturnCodes::EVENT_OK;
0318 }
0319 
0320 void RawClusterBuilderTemplate::PrintCylGeom(RawTowerGeomContainer *towergeom, const std::string &fname) const
0321 {
0322   std::ofstream outfile(fname);
0323   if (!outfile.is_open())
0324   {
0325     std::cout << "Error in BEmcRec::RawClusterBuilderTemplate::PrintCylGeom(): Failed to open file "
0326               << fname << std::endl;
0327     return;
0328   }
0329   outfile << NBINX << " " << NBINY << std::endl;
0330   for (int ip = 0; ip < NBINX; ip++)
0331   {
0332     outfile << ip << " " << towergeom->get_phicenter(ip) << std::endl;
0333   }
0334   for (int ip = 0; ip < NBINY; ip++)
0335   {
0336     outfile << ip << " " << towergeom->get_etacenter(ip) << std::endl;
0337   }
0338   outfile.close();
0339 }
0340 
0341 bool RawClusterBuilderTemplate::Cell2Abs(RawTowerGeomContainer * /*towergeom*/, float /*phiC*/, float /*etaC*/, float &phi, float &eta)
0342 {
0343   phi = eta = 0;
0344   return false;
0345 }
0346 
0347 int RawClusterBuilderTemplate::process_event(PHCompositeNode *topNode)
0348 {
0349   if (bemc == nullptr)
0350   {
0351     std::cout << "Error in RawClusterBuilderTemplate::process_event(): detector is not defined; use RawClusterBuilderTemplate::Detector() to define it" << std::endl;
0352     return Fun4AllReturnCodes::ABORTEVENT;
0353   }
0354 
0355   RawTowerContainer *towers = nullptr;
0356   if (m_UseTowerInfo < 1)
0357   {
0358     std::string towernodename = "TOWER_CALIB_" + detector;
0359 
0360     // Grab the towers
0361     towers = findNode::getClass<RawTowerContainer>(topNode, towernodename);
0362     if (!towers)
0363     {
0364       std::cout << PHWHERE << ": Could not find node " << towernodename << std::endl;
0365       return Fun4AllReturnCodes::DISCARDEVENT;
0366     }
0367   }
0368 
0369   // At this stage, it is more efficient to read the simple geometry node in any case
0370   // Indeed, we only need to read the calorimeter ID
0371   m_TowerGeomNodeName = "TOWERGEOM_" + detector;
0372   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, m_TowerGeomNodeName);
0373   if (!towergeom)
0374   {
0375     std::cout << PHWHERE << ": Could not find node " << m_TowerGeomNodeName << std::endl;
0376     return Fun4AllReturnCodes::ABORTEVENT;
0377   }
0378   TowerInfoContainer *calib_towerinfos = nullptr;
0379   if (m_UseTowerInfo > 0)
0380   {
0381     std::string towerinfoNodename = "TOWERINFO_CALIB_" + detector;
0382     if (!m_inputnodename.empty())
0383     {
0384       towerinfoNodename = m_inputnodename;
0385     }
0386 
0387     calib_towerinfos = findNode::getClass<TowerInfoContainer>(topNode, towerinfoNodename);
0388     if (!calib_towerinfos)
0389     {
0390       std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
0391                 << " " << towerinfoNodename << " Node missing, doing bail out!"
0392                 << std::endl;
0393 
0394       return Fun4AllReturnCodes::DISCARDEVENT;
0395     }
0396   }
0397 
0398   // Get vertex
0399   float vx = 0;
0400   float vy = 0;
0401   float vz = 0;
0402   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0403 
0404   if (vertexmap && m_UseAltZVertex == 0)  // default
0405   {
0406     GlobalVertex* vtx = vertexmap->begin()->second;
0407     if (vtx)
0408     {
0409       auto typeStartIter = vtx->find_vertexes(m_vertex_type);
0410       auto typeEndIter = vtx->end_vertexes();
0411       for (auto iter = typeStartIter; iter != typeEndIter; ++iter)
0412       {
0413         const auto& [type, vertexVec] = *iter;
0414         if (type != m_vertex_type)
0415         {
0416           continue;
0417         }
0418         for (const auto* vertex : vertexVec)
0419         {
0420           if (!vertex)
0421           {
0422             continue;
0423           }
0424           vx = vertex->get_x();
0425           vy = vertex->get_y();
0426           vz = vertex->get_z();
0427         }
0428       }
0429     }
0430   }
0431 
0432   MbdVertexMap *mbdmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0433 
0434   if (mbdmap && m_UseAltZVertex == 1)
0435   {
0436     MbdVertex *bvertex = nullptr;
0437     for (MbdVertexMap::ConstIter mbditer = mbdmap->begin();
0438          mbditer != mbdmap->end();
0439          ++mbditer)
0440     {
0441       bvertex = mbditer->second;
0442     }
0443     if (bvertex)
0444     {
0445       vz = bvertex->get_z();
0446     }
0447   }
0448   
0449   if (m_UseAltZVertex == 3)
0450   {
0451     PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0452     if (truthinfo)
0453     {
0454       PHG4TruthInfoContainer::VtxRange vtxrange = truthinfo->GetVtxRange();
0455       for (PHG4TruthInfoContainer::ConstVtxIterator iter = vtxrange.first; iter != vtxrange.second; ++iter) 
0456       {
0457          PHG4VtxPoint *vtx_tr = iter->second;
0458          if ( vtx_tr->get_id() == 1 )
0459          { 
0460            vz = vtx_tr->get_z();
0461            vy = vtx_tr->get_y();
0462            vx = vtx_tr->get_x();
0463          }
0464       }
0465     }
0466     else {
0467       std::cout << "RawClusterBuilderTemplate: Error requiring truth vertex but non was found. Exiting" << std::endl;  
0468       return Fun4AllReturnCodes::ABORTEVENT;
0469     } 
0470   }
0471 
0472   // Set vertex
0473   float vertex[3] = {vx, vy, vz};
0474   bemc->SetVertex(vertex);
0475   // Set threshold
0476   bemc->SetTowerThreshold(_min_tower_e);
0477   bemc->SetPeakThreshold(_min_peak_e);
0478 
0479   bemc->SetProbNoiseParam(fProbNoiseParam);
0480   bemc->SetProfileProb(bProfProb);
0481 
0482   _clusters->Reset();  // make sure cluster container is empty before filling it with new clusters
0483 
0484   // Define vector of towers in EmcModule format to input into BEmc
0485   EmcModule vhit;
0486   std::vector<EmcModule> HitList;
0487   HitList.erase(HitList.begin(), HitList.end());
0488   int ich;
0489 
0490   if (m_UseTowerInfo < 1)
0491   {
0492     // make the list of towers above threshold
0493     RawTowerContainer::ConstRange begin_end = towers->getTowers();
0494     RawTowerContainer::ConstIterator itr = begin_end.first;
0495 
0496     for (; itr != begin_end.second; ++itr)
0497     {
0498       RawTower *tower = itr->second;
0499       // debug change
0500       //      std::cout << "  Tower e = " << tower->get_energy()
0501       //           << " (" << _min_tower_e << ")" << std::endl;
0502       if (IsAcceptableTower(tower))
0503       {
0504         // debug change
0505         // std::cout << "(" << tower->get_column() << "," << tower->get_row()
0506         //      << ")  (" << tower->get_binphi() << "," << tower->get_bineta()
0507         //      << ")" << std::endl;
0508         //    ix = tower->get_column();
0509         RawTowerDefs::keytype towerid = tower->get_id();
0510         int ix = RawTowerDefs::decode_index2(towerid);  // index2 is phi in CYL
0511         int iy = RawTowerDefs::decode_index1(towerid);  // index1 is eta in CYL
0512         ix -= BINX0;
0513         iy -= BINY0;
0514         // debug change
0515         //      ix = tower->get_bineta() - BINX0;  // eta: index1
0516         //      iy = tower->get_binphi() - BINY0;  // phi: index2
0517         if (ix >= 0 && ix < NBINX && iy >= 0 && iy < NBINY)
0518         {
0519           ich = iy * NBINX + ix;
0520           vhit.ich = ich;
0521           vhit.amp = tower->get_energy() * fEnergyNorm;  // !!! Global Calibration
0522           vhit.tof = 0.;
0523           HitList.push_back(vhit);
0524         }
0525       }
0526     }
0527   }
0528   else if (m_UseTowerInfo)
0529   {
0530     // make the list of towers above threshold
0531     // debug change
0532     // TowerInfoContainer::ConstRange begin_end = calib_towerinfos->getTowers();
0533     // TowerInfoContainer::ConstIterator rtiter;
0534     unsigned int nchannels = calib_towerinfos->size();
0535     // for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0536 
0537     //float total_energy = 0;
0538     for (unsigned int channel = 0; channel < nchannels; channel++)
0539     {
0540       TowerInfo *tower_info = calib_towerinfos->get_tower_at_channel(channel);
0541       
0542       // debug change
0543       //      std::cout << "  Tower e = " << tower->get_energy()
0544       //           << " (" << _min_tower_e << ")" << std::endl;
0545       if (IsAcceptableTower(tower_info))
0546       {
0547         unsigned int towerkey = calib_towerinfos->encode_key(channel);
0548         int ieta = calib_towerinfos->getTowerEtaBin(towerkey);
0549         int iphi = calib_towerinfos->getTowerPhiBin(towerkey);
0550 
0551         const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0552 
0553         int ix = RawTowerDefs::decode_index2(key);  // index2 is phi in CYL
0554         int iy = RawTowerDefs::decode_index1(key);  // index1 is eta in CYL
0555 
0556         ix -= BINX0;
0557         iy -= BINY0;
0558 
0559         if (ix >= 0 && ix < NBINX && iy >= 0 && iy < NBINY)
0560         {
0561           ich = iy * NBINX + ix;
0562           // add key field to vhit
0563           vhit.ich = ich;
0564           vhit.amp = tower_info->get_energy() * fEnergyNorm;  // !!! Global Calibration
0565           vhit.tof = tower_info->get_time();
0566 
0567           // debug change
0568           /*total_energy += vhit.amp;
0569           
0570           if (vhit.amp > 1e-8) {
0571             std::cout << "hit: (" << iy << ", " << ix << ", " << vhit.amp << ", " << tower_info->get_energy() << ", " << vhit.tof << ")\n";
0572           }
0573           */
0574 
0575           HitList.push_back(vhit);
0576         }
0577       }
0578     }
0579     // debug change
0580     //std::cout << "Total hit energy = " << total_energy << std::endl;
0581   }
0582 
0583   bemc->SetModules(&HitList);
0584 
0585   // Find clusters (as a set of towers with common edge)
0586   int ncl = bemc->FindClusters();
0587   if (ncl < 0)
0588   {
0589     std::cout << "!!! Error in BEmcRec::FindClusters(): numbers of cluster "
0590               << ncl << " ?" << std::endl;
0591     return Fun4AllReturnCodes::ABORTEVENT;
0592   }
0593 
0594   // Get pointer to clusters
0595   std::vector<EmcCluster> *ClusterList = bemc->GetClusters();
0596   std::vector<EmcCluster>::iterator pc;
0597 
0598   std::vector<EmcCluster>::iterator pp;
0599   float ecl;
0600   float ecore;
0601   float xcg;
0602   float ycg;
0603   float xx;
0604   float xy;
0605   float yy;
0606   //  float xcorr, ycorr;
0607   EmcModule hmax;
0608   RawCluster *cluster;
0609 
0610   std::vector<EmcCluster> PList;
0611   std::vector<EmcModule> Peaks;
0612 
0613   float prob;
0614   float chi2;
0615   int ndf;
0616   float xg;
0617   float yg;
0618   float zg;
0619 
0620   std::vector<EmcModule>::iterator ph;
0621   std::vector<EmcModule> hlist;
0622 
0623   // ncl = 0;
0624   for (pc = ClusterList->begin(); pc != ClusterList->end(); ++pc)
0625   {
0626     //    ecl = pc->GetTotalEnergy();
0627     //    pc->GetMoments( &xcg, &ycg, &xx, &xy, &yy );
0628 
0629     int npk = pc->GetSubClusters(PList, Peaks,m_subclustersplitting);
0630     if (npk < 0)
0631     {
0632       return Fun4AllReturnCodes::ABORTEVENT;
0633     }
0634 
0635     // debug change
0636     //    std::cout << "  iCl = " << ncl << " (" << npk << "): E ="
0637     //         << ecl << "  x = " << xcg << "  y = " << ycg << std::endl;
0638 
0639     for (pp = PList.begin(); pp != PList.end(); ++pp)
0640     {
0641       // Cluster energy
0642       ecl = pp->GetTotalEnergy();
0643       if (ecl < m_min_cluster_e)
0644       {
0645         continue;
0646       }
0647       ecore = pp->GetECoreCorrected();
0648       // 3x3 energy around center of gravity
0649       // e9 = pp->GetE9();
0650       // Ecore (basically near 2x2 energy around center of gravity)
0651       // ecore = pp->GetECore();
0652       // Center of Gravity etc.
0653       pp->GetMoments(xcg, ycg, xx, xy, yy);
0654 
0655       if (m_UseAltZVertex == 2)
0656       {
0657         xg = -99999;  // signal to force zvtx = 0
0658         pp->GetGlobalPos(xg, yg, zg);
0659       }
0660       else
0661       {
0662         xg = 0;  // usual mode, uses regular zvtx
0663         pp->GetGlobalPos(xg, yg, zg);
0664       }
0665 
0666       // Tower with max energy
0667       hmax = pp->GetMaxTower();
0668 
0669       // debug change
0670       //      phi = (xcg-float(NPHI)/2.+0.5)/float(NPHI)*2.*M_PI;
0671       //      eta = (ycg-float(NETA)/2.+0.5)/float(NETA)*2.2; // -1.1<eta<1.1;
0672 
0673       //      Cell2Abs(towergeom,xcg,ycg,phi,eta);
0674 
0675       //      pp->GetCorrPos(&xcorr, &ycorr);
0676       //      Cell2Abs(towergeom, xcorr, ycorr, phi, eta);
0677       //      const double ref_radius = towergeom->get_radius();
0678 
0679       //      phi = 0;
0680       //      if (phi > M_PI) phi -= 2. * M_PI;  // convert to [-pi,pi]]
0681 
0682       //      prob = -1;
0683       chi2 = 0;
0684       ndf = 0;
0685       prob = pp->GetProb(chi2, ndf);
0686       //      std::cout << "Prob/Chi2/NDF = " << prob << " " << chi2
0687       //           << " " << ndf << " Ecl = " << ecl << std::endl;
0688 
0689       cluster = new RawClusterv1();
0690       cluster->set_energy(ecl);
0691       cluster->set_ecore(ecore);
0692       cluster->set_r(std::sqrt(xg * xg + yg * yg));
0693       cluster->set_phi(std::atan2(yg, xg));
0694       cluster->set_z(zg);
0695 
0696       cluster->set_prob(prob);
0697       if (ndf > 0)
0698       {
0699         cluster->set_chi2(chi2 / ndf);
0700       }
0701       else
0702       {
0703         cluster->set_chi2(0);
0704       }
0705       hlist = pp->GetHitList();
0706       ph = hlist.begin();
0707 
0708       // accumulate energy-weighted time
0709       float ew_num = 0.0;   // sum(E * t)
0710       float ew_den = 0.0;   // sum(E)
0711       bool   saw_nonzero_t = false;
0712 
0713       while (ph != hlist.end())
0714       {
0715         ich = (*ph).ich;
0716         int iy = ich / NBINX;
0717         int ix = ich % NBINX;
0718         // that code needs a closer look - here are the towers
0719         // with their energy added to the cluster object where
0720         // the id is the tower id
0721         // !!!!! Make sure twrkey is correctly extracted
0722         //        RawTowerDefs::keytype twrkey = RawTowerDefs::encode_towerid(towers->getCalorimeterID(), ix + BINX0, iy + BINY0);
0723         const RawTowerDefs::CalorimeterId Calo_ID = towergeom->get_calorimeter_id();
0724         RawTowerDefs::keytype twrkey = RawTowerDefs::encode_towerid(Calo_ID, iy + BINY0, ix + BINX0);  // Becuase in this part index1 is iy
0725         //  std::cout << iphi << " " << ieta << ": "
0726         //           << twrkey << " e = " << (*ph).amp) << std::endl;
0727           
0728         const float amp = (*ph).amp;
0729         const float tof = (*ph).tof;
0730 
0731         // add tower (energy is un-normalized inside vhit)
0732         cluster->addTower(twrkey, amp / fEnergyNorm);
0733 
0734         // accumulate EW time (finite guard; treat exact 0 as “no time”)
0735         if (std::isfinite(tof))
0736         {
0737           if (std::abs(tof) > 1e-9F) { saw_nonzero_t = true; }
0738           ew_num += amp * tof;    // float math end-to-end
0739         }
0740         ew_den += amp;
0741 
0742 
0743         ++ph;
0744       }
0745 
0746       // stamp tower CoG (raw & corrected) into the cluster (now via properties in v1)
0747       float xcorr = xcg;
0748       float ycorr = ycg;
0749       bemc->CorrectPosition(ecl, xcg, ycg, xcorr, ycorr);
0750       cluster->set_tower_cog(xcg, ycg, xcorr, ycorr);
0751 
0752       const float tmean = (ew_den > 0.0F && saw_nonzero_t)
0753                             ? (ew_num / ew_den)
0754                             : std::numeric_limits<float>::quiet_NaN();
0755       cluster->set_mean_time(tmean);
0756         // mechanical incidence angles (CEMC-only; stored as signed alphas in properties)
0757         {
0758           BEmcRecCEMC* cemc = dynamic_cast<BEmcRecCEMC*>(bemc);
0759           if (cemc)
0760           {
0761             float a_phi_sgn = std::numeric_limits<float>::quiet_NaN();
0762             float a_eta_sgn = std::numeric_limits<float>::quiet_NaN();
0763 
0764             // use corrected CoG for incidence, in tower units
0765             const bool inc_ok = cemc->calculateIncidence(
0766                 xcorr, ycorr,
0767                 a_phi_sgn, a_eta_sgn);
0768 
0769             if (inc_ok)
0770             {
0771               // Use explicit property IDs so this compiles even when the
0772               // RawCluster header does not yet define the named incidence enums.
0773               cluster->set_property(RawCluster::prop_incidence_alpha_phi, a_phi_sgn);
0774               cluster->set_property(RawCluster::prop_incidence_alpha_eta, a_eta_sgn);
0775             }
0776           }
0777         }
0778 
0779       _clusters->AddCluster(cluster);
0780       // ncl++;
0781 
0782       //      std::cout << "    ipk = " << ipk << ": E = " << ecl << "  E9 = "
0783       //           << e9 << "  x = " << xcg << "  y = " << ycg
0784       //           << "  MaxTower: (" << hmax.ich%NPHI << ","
0785       //           << hmax.ich/NPHI << ") e = " << hmax.amp << std::endl;
0786     }
0787   }
0788 
0789   if (chkenergyconservation && towers && _clusters)
0790   {
0791     double ecluster = _clusters->getTotalEdep();
0792     double etower = towers->getTotalEdep();
0793     if (ecluster > 0)
0794     {
0795       if (fabs(etower - ecluster) / ecluster > 1e-9)
0796       {
0797         std::cout << "energy conservation violation: ETower: " << etower
0798                   << " ECluster: " << ecluster
0799                   << " diff: " << etower - ecluster << std::endl;
0800       }
0801     }
0802     else
0803     {
0804       if (etower != 0)
0805       {
0806         std::cout << "energy conservation violation: ETower: " << etower
0807                   << " ECluster: " << ecluster << std::endl;
0808       }
0809     }
0810   }
0811   else if (chkenergyconservation)
0812   {
0813     std::cout << "RawClusterBuilderTemplate : energy conservation check asked for but tower or cluster container is NULL" << std::endl;
0814   }
0815 
0816   return Fun4AllReturnCodes::EVENT_OK;
0817 }
0818 
0819 void RawClusterBuilderTemplate::CreateNodes(PHCompositeNode *topNode)
0820 {
0821   PHNodeIterator iter(topNode);
0822 
0823   // Grab the cEMC node
0824   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0825   if (!dstNode)
0826   {
0827     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0828     throw std::runtime_error("Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0829   }
0830 
0831   // Get the _det_name subnode
0832   PHCompositeNode *cemcNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", detector));
0833 
0834   // Check that it is there
0835   if (!cemcNode)
0836   {
0837     cemcNode = new PHCompositeNode(detector);
0838     dstNode->addNode(cemcNode);
0839   }
0840   ClusterNodeName = "CLUSTER_" + detector;
0841   if (!m_outputnodename.empty())
0842   {
0843     ClusterNodeName = m_outputnodename;
0844   }
0845   else if (m_UseTowerInfo)
0846   {
0847     ClusterNodeName = "CLUSTERINFO_" + detector;
0848   }
0849   _clusters = findNode::getClass<RawClusterContainer>(dstNode, ClusterNodeName);
0850   if (!_clusters)
0851   {
0852     _clusters = new RawClusterContainer();
0853   }
0854 
0855   PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_clusters, ClusterNodeName, "PHObject");
0856   cemcNode->addNode(clusterNode);
0857 }
0858 
0859 bool RawClusterBuilderTemplate::IsAcceptableTower(TowerInfo *tower) const
0860 {
0861   if (tower->get_energy() < _min_tower_e)
0862   {
0863     return false;
0864   }
0865 
0866   if (m_do_tower_selection)
0867   {
0868     if (!tower->get_isGood())
0869     {
0870       return false;
0871     }
0872 
0873   }
0874   return true;
0875 }
0876 
0877 bool RawClusterBuilderTemplate::IsAcceptableTower(RawTower *tower) const
0878 {
0879   if (tower->get_energy() < _min_tower_e)
0880   {
0881     return false;
0882   }
0883   return true;
0884 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!