Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:17:33 

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