packages/trigger/CaloTriggerSim.cc

0001 #include "CaloTriggerSim.h"
0002
0003 #include "CaloTriggerInfo.h"
0004 #include "CaloTriggerInfov1.h"
0005
0006 // sPHENIX includes
0007 #include <calobase/RawTowerDefs.h>  // for encode_towerid
0008 #include <calobase/RawTowerGeom.h>
0009 #include <calobase/RawTowerGeomContainer.h>
0010 #include <calobase/TowerInfo.h>
0011 #include <calobase/TowerInfoContainer.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/SubsysReco.h>
0015
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHIODataNode.h>
0018 #include <phool/PHNode.h>
0019 #include <phool/PHNodeIterator.h>
0020 #include <phool/PHObject.h>
0021 #include <phool/getClass.h>
0022 #include <phool/phool.h>
0023
0024 // standard includes
0025 #include <algorithm>
0026 #include <cmath>
0027 #include <cstdlib>
0028 #include <iomanip>
0029 #include <iostream>
0030 #include <utility>
0031 #include <vector>
0032
0033 CaloTriggerSim::CaloTriggerSim(const std::string &name)
0034   : SubsysReco(name)
0035 {
0036   return;
0037 }
0038
0039 double CaloTriggerSim::truncate_8bit(const double raw_E)
0040 {
0041   double rawE = std::min(raw_E, 45.0);
0042   int counts = std::floor(rawE / (45.0 / 256));
0043
0044   return counts * (45.0 / 256);
0045 }
0046 int CaloTriggerSim::InitRun(PHCompositeNode *topNode)
0047 {
0048   return CreateNode(topNode);
0049 }
0050
0051 int CaloTriggerSim::process_event(PHCompositeNode *topNode)
0052 {
0053   if (Verbosity() > 0)
0054   {
0055     std::cout << "CaloTriggerSim::process_event: entering" << std::endl;
0056   }
0057
0058   // pull out the tower containers and geometry objects at the start
0059
0060   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0061   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0062   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0063
0064   if (Verbosity() > 0)
0065   {
0066     std::cout << "CaloTriggerSim::process_event: " << towersEM3->size() << " TOWERINFO_CALIB_CEMC towers" << std::endl;
0067     std::cout << "CaloTriggerSim::process_event: " << towersIH3->size() << " TOWERINFO_CALIB_HCALIN towers" << std::endl;
0068     std::cout << "CaloTriggerSim::process_event: " << towersOH3->size() << " TOWERINFO_CALIB_HCALOUT towers" << std::endl;
0069   }
0070
0071   RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0072   RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0073   RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0074
0075   // get the binning from the geometry (different for 1D vs 2D...)
0076   int geom_etabins = geomEM->get_etabins();
0077   int geom_phibins = geomEM->get_phibins();
0078
0079   // if internal knowledge of geometry is unset, set it now (should
0080   // only happen once, on the first event)
0081   if (m_EMCAL_1x1_NETA < 0)
0082   {
0083     m_EMCAL_1x1_NETA = geom_etabins;
0084     m_EMCAL_1x1_NPHI = geom_phibins;
0085
0086     // half as many 2x2 windows along each axis as 1x1
0087     m_EMCAL_2x2_NETA = geom_etabins / 2;
0088     m_EMCAL_2x2_NPHI = geom_phibins / 2;
0089
0090     // each 2x2 window defines a 4x4 window for which that 2x2 window
0091     // is the upper-left corner, so there are as many 4x4's as 2x2's
0092     // (except in eta, where the edge effect means there is 1 fewer)
0093     m_EMCAL_4x4_NETA = geom_etabins / 2 - 1;
0094     m_EMCAL_4x4_NPHI = geom_phibins / 2;
0095
0096     // reset all maps
0097     m_EMCAL_1x1_MAP.resize(m_EMCAL_1x1_NETA, std::vector<double>(m_EMCAL_1x1_NPHI, 0));
0098     m_EMCAL_2x2_MAP.resize(m_EMCAL_2x2_NETA, std::vector<double>(m_EMCAL_2x2_NPHI, 0));
0099     m_EMCAL_4x4_MAP.resize(m_EMCAL_4x4_NETA, std::vector<double>(m_EMCAL_4x4_NPHI, 0));
0100
0101     if (Verbosity() > 0)
0102     {
0103       std::cout << "CaloTriggerSim::process_event: setting number of window in eta / phi,";
0104       std::cout << "1x1 are " << m_EMCAL_1x1_NETA << " / " << m_EMCAL_1x1_NPHI << ", ";
0105       std::cout << "2x2 are " << m_EMCAL_2x2_NETA << " / " << m_EMCAL_2x2_NPHI << ", ";
0106       std::cout << "4x4 are " << m_EMCAL_4x4_NETA << " / " << m_EMCAL_4x4_NPHI << std::endl;
0107     }
0108   }
0109
0110   // reset 1x1 map
0111   fill(m_EMCAL_1x1_MAP.begin(), m_EMCAL_1x1_MAP.end(), std::vector<double>(m_EMCAL_1x1_NPHI, 0));
0112
0113   // iterate over EMCal towers, constructing 1x1's
0114   unsigned int ntowers_EM = towersEM3->size();
0115   for (unsigned int channel = 0; channel < ntowers_EM; channel++)
0116   {
0117     TowerInfo *tower = towersEM3->get_tower_at_channel(channel);
0118
0119     // RawTowerGeom *tower_geom = geomEM->get_tower_geometry(tower->get_key());
0120
0121     unsigned int towerkey = towersEM3->encode_key(channel);
0122     int this_etabin = towersEM3->getTowerEtaBin(towerkey);
0123     int this_phibin = towersEM3->getTowerPhiBin(towerkey);
0124     double this_E = tower->get_energy();
0125
0126     m_EMCAL_1x1_MAP[this_etabin][this_phibin] += this_E;
0127
0128     if (Verbosity() > 1 && tower->get_energy() > 1)
0129     {
0130       std::cout << "CaloTriggerSim::process_event: EMCal 1x1 tower eta  bin  / phi  bin  / E = " << std::setprecision(6) << this_etabin << "  / " << this_phibin << " / " << this_E << std::endl;
0131     }
0132   }
0133
0134   // reset 2x2 map and best
0135   fill(m_EMCAL_2x2_MAP.begin(), m_EMCAL_2x2_MAP.end(), std::vector<double>(m_EMCAL_2x2_NPHI, 0));
0136
0137   m_EMCAL_2x2_BEST_E = 0;
0138   m_EMCAL_2x2_BEST_PHI = 0;
0139   m_EMCAL_2x2_BEST_ETA = 0;
0140
0141   // now reconstruct 2x2 map from 1x1 map
0142   for (std::vector<double>::size_type ieta = 0; ieta < m_EMCAL_2x2_NETA; ieta++)
0143   {
0144     for (std::vector<double>::size_type iphi = 0; iphi < m_EMCAL_2x2_NPHI; iphi++)
0145     {
0146       double this_sum = 0;
0147
0148       this_sum += m_EMCAL_1x1_MAP[2 * ieta][2 * iphi];
0149       this_sum += m_EMCAL_1x1_MAP[2 * ieta][(2 * iphi) + 1];  // 2 * iphi + 1 is safe, since m_EMCAL_2x2_NPHI = m_EMCAL_1x1_NPHI / 2
0150       this_sum += m_EMCAL_1x1_MAP[(2 * ieta) + 1][2 * iphi];  // 2 * ieta + 1 is safe, since m_EMCAL_2x2_NETA = m_EMCAL_1x1_NETA / 2
0151       this_sum += m_EMCAL_1x1_MAP[(2 * ieta) + 1][(2 * iphi) + 1];
0152
0153       if (m_EmulateTruncationFlag)
0154       {
0155         this_sum = truncate_8bit(this_sum);
0156       }
0157
0158       // populate 2x2 map
0159       m_EMCAL_2x2_MAP[ieta][iphi] = this_sum;
0160
0161       // to calculate the eta, phi position, take the average of that of the 1x1's
0162       double this_eta = 0.5 * (geomEM->get_etacenter(2 * ieta) + geomEM->get_etacenter((2 * ieta) + 1));
0163       double this_phi = 0.5 * (geomEM->get_phicenter(2 * iphi) + geomEM->get_phicenter((2 * iphi) + 1));
0164       // wrap-around phi (apparently needed for 2D geometry?)
0165       if (this_phi > M_PI)
0166       {
0167         this_phi -= 2 * M_PI;
0168       }
0169       if (this_phi < -M_PI)
0170       {
0171         this_phi += 2 * M_PI;
0172       }
0173
0174       if (Verbosity() > 1 && this_sum > 1)
0175       {
0176         std::cout << "CaloTriggerSim::process_event: EMCal 2x2 tower eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << ieta << " ) / " << this_phi << " ( " << iphi << " ) / " << this_sum << std::endl;
0177       }
0178
0179       if (this_sum > m_EMCAL_2x2_BEST_E)
0180       {
0181         m_EMCAL_2x2_BEST_E = this_sum;
0182         m_EMCAL_2x2_BEST_PHI = this_phi;
0183         m_EMCAL_2x2_BEST_ETA = this_eta;
0184       }
0185     }
0186   }
0187
0188   if (Verbosity() > 0)
0189   {
0190     std::cout << "CaloTriggerSim::process_event: best EMCal 2x2 window is at eta / phi = " << m_EMCAL_2x2_BEST_ETA << " / " << m_EMCAL_2x2_BEST_PHI << " and E = " << m_EMCAL_2x2_BEST_E << std::endl;
0191   }
0192
0193   // reset 4x4 map & best
0194   fill(m_EMCAL_4x4_MAP.begin(), m_EMCAL_4x4_MAP.end(), std::vector<double>(m_EMCAL_4x4_NPHI, 0));
0195
0196   m_EMCAL_4x4_BEST_E = 0;
0197   m_EMCAL_4x4_BEST_PHI = 0;
0198   m_EMCAL_4x4_BEST_ETA = 0;
0199
0200   int emcal_4x4_best_iphi = -1;
0201   int emcal_4x4_best_ieta = -1;
0202
0203   // now reconstruct (sliding) 4x4 map from 2x2 map
0204   for (int ieta = 0; ieta < m_EMCAL_4x4_NETA; ieta++)
0205   {
0206     for (int iphi = 0; iphi < m_EMCAL_4x4_NPHI; iphi++)
0207     {
0208       // for eta calculation (since eta distribution is potentially
0209       // non-uniform), average positions of all four towers
0210       double this_eta = 0.25 * (geomEM->get_etacenter(2 * ieta) + geomEM->get_etacenter((2 * ieta) + 1) + geomEM->get_etacenter((2 * ieta) + 2) + geomEM->get_etacenter((2 * ieta) + 3));
0211       // for phi calculation (since phi distribution is uniform), take
0212       // first tower and add 1.5 tower widths
0213       double this_phi = geomEM->get_phicenter(2 * iphi) + (1.5 * (geomEM->get_phicenter((2 * iphi) + 1) - geomEM->get_phicenter(2 * iphi)));
0214       // wrap-around phi (apparently needed for 2D geometry?)
0215       if (this_phi > M_PI)
0216       {
0217         this_phi -= 2 * M_PI;
0218       }
0219       if (this_phi < -M_PI)
0220       {
0221         this_phi += 2 * M_PI;
0222       }
0223
0224       double this_sum = 0;
0225
0226       this_sum += m_EMCAL_2x2_MAP[ieta][iphi];
0227       this_sum += m_EMCAL_2x2_MAP[ieta + 1][iphi];  // ieta + 1 is safe, since m_EMCAL_4x4_NETA = m_EMCAL_2x2_NETA - 1
0228
0229       if (iphi != m_EMCAL_4x4_NPHI - 1)
0230       {
0231         // if we are not in the last phi row, can safely access 'iphi+1'
0232         this_sum += m_EMCAL_2x2_MAP[ieta][iphi + 1];
0233         this_sum += m_EMCAL_2x2_MAP[ieta + 1][iphi + 1];
0234       }
0235       else
0236       {
0237         // if we are in the last phi row, wrap back around to zero
0238         this_sum += m_EMCAL_2x2_MAP[ieta][0];
0239         this_sum += m_EMCAL_2x2_MAP[ieta + 1][0];
0240       }
0241
0242       m_EMCAL_4x4_MAP[ieta][iphi] = this_sum;
0243
0244       if (Verbosity() > 1 && this_sum > 1)
0245       {
0246         std::cout << "CaloTriggerSim::process_event: EMCal 4x4 tower eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << ieta << " ) / " << this_phi << " ( " << iphi << " ) / " << this_sum << std::endl;
0247       }
0248
0249       if (this_sum > m_EMCAL_4x4_BEST_E)
0250       {
0251         m_EMCAL_4x4_BEST_E = this_sum;
0252         m_EMCAL_4x4_BEST_PHI = this_phi;
0253         m_EMCAL_4x4_BEST_ETA = this_eta;
0254
0255         emcal_4x4_best_iphi = iphi;
0256         emcal_4x4_best_ieta = ieta;
0257       }
0258     }
0259   }
0260
0261   m_EMCAL_4x4_BEST2_E = 0;
0262   m_EMCAL_4x4_BEST2_PHI = 0;
0263   m_EMCAL_4x4_BEST2_ETA = 0;
0264
0265   // find second-largest 4x4 which is > 1 tower away...
0266   for (int ieta = 0; ieta < m_EMCAL_4x4_NETA; ieta++)
0267   {
0268     for (int iphi = 0; iphi < m_EMCAL_4x4_NPHI; iphi++)
0269     {
0270       int deta = ieta - emcal_4x4_best_ieta;
0271       int dphi = (iphi - emcal_4x4_best_iphi) % m_EMCAL_4x4_NPHI;
0272
0273       if (abs(deta) < 1.5 && abs(dphi) < 1.5)
0274       {
0275         continue;
0276       }
0277
0278       double this_eta = 0.25 * (geomEM->get_etacenter(2 * ieta) + geomEM->get_etacenter((2 * ieta) + 1) + geomEM->get_etacenter((2 * ieta) + 2) + geomEM->get_etacenter((2 * ieta) + 3));
0279       double this_phi = geomEM->get_phicenter(2 * iphi) + (1.5 * (geomEM->get_phicenter((2 * iphi) + 1) - geomEM->get_phicenter(2 * iphi)));
0280
0281       if (this_phi > M_PI)
0282       {
0283         this_phi -= 2 * M_PI;
0284       }
0285       if (this_phi < -M_PI)
0286       {
0287         this_phi += 2 * M_PI;
0288       }
0289
0290       double this_sum = m_EMCAL_4x4_MAP[ieta][iphi];
0291
0292       if (this_sum > m_EMCAL_4x4_BEST2_E)
0293       {
0294         m_EMCAL_4x4_BEST2_E = this_sum;
0295         m_EMCAL_4x4_BEST2_PHI = this_phi;
0296         m_EMCAL_4x4_BEST2_ETA = this_eta;
0297       }
0298     }
0299   }
0300
0301   if (Verbosity() > 0)
0302   {
0303     std::cout << "CaloTriggerSim::process_event: best EMCal 4x4 window is at eta / phi = " << m_EMCAL_4x4_BEST_ETA << " / " << m_EMCAL_4x4_BEST_PHI << " and E = " << m_EMCAL_4x4_BEST_E << std::endl;
0304     std::cout << "CaloTriggerSim::process_event: 2nd best EMCal 4x4 window is at eta / phi = " << m_EMCAL_4x4_BEST2_ETA << " / " << m_EMCAL_4x4_BEST2_PHI << " and E = " << m_EMCAL_4x4_BEST2_E << std::endl;
0305   }
0306
0307   // begin full calo sim
0308
0309   // get the 0.1x0.1 binning from the OHCal geometry
0310   int geomOH_etabins = geomOH->get_etabins();
0311   int geomOH_phibins = geomOH->get_phibins();
0312
0313   // if internal knowledge of geometry is unset, set it now
0314   if (m_FULLCALO_0p1x0p1_NETA < 0)
0315   {
0316     m_FULLCALO_PHI_START = geomOH->get_phibounds(0).first;
0317     m_FULLCALO_PHI_END = geomOH->get_phibounds(geomOH_phibins - 1).second;
0318
0319     m_FULLCALO_0p1x0p1_NETA = geomOH_etabins;
0320     m_FULLCALO_0p1x0p1_NPHI = geomOH_phibins;
0321
0322     // half as many 0.2x0.2 windows along each axis as 0.1x0.1
0323     m_FULLCALO_0p2x0p2_NETA = geomOH_etabins / 2;
0324     m_FULLCALO_0p2x0p2_NPHI = geomOH_phibins / 2;
0325
0326     // each 0.2x0.2 window defines a 0.4x0.4 window for which that
0327     // 0.2x0.2 window is the upper-left corner, so there are as many
0328     // 0.4x0.4's as 0.2x0.2's (except in eta, where the edge effect
0329     // means there is 1 fewer)
0330     m_FULLCALO_0p4x0p4_NETA = geomOH_etabins / 2 - 1;
0331     m_FULLCALO_0p4x0p4_NPHI = geomOH_phibins / 2;
0332
0333     // for 0.6x0.6 windows, the above logic applies, except that the
0334     // edge effect causes there to be 2 fewer less in eta
0335     m_FULLCALO_0p6x0p6_NETA = geomOH_etabins / 2 - 2;
0336     m_FULLCALO_0p6x0p6_NPHI = geomOH_phibins / 2;
0337
0338     // for 0.8x0.8 windows, the above logic applies, except that the
0339     // edge effect causes there to be 3 fewer less in eta
0340     m_FULLCALO_0p8x0p8_NETA = geomOH_etabins / 2 - 3;
0341     m_FULLCALO_0p8x0p8_NPHI = geomOH_phibins / 2;
0342
0343     // for 1.0x1.0 windows, the above logic applies, except that the
0344     // edge effect causes there to be 4 fewer less in eta
0345     m_FULLCALO_1p0x1p0_NETA = geomOH_etabins / 2 - 4;
0346     m_FULLCALO_1p0x1p0_NPHI = geomOH_phibins / 2;
0347
0348     // reset all maps
0349     m_FULLCALO_0p1x0p1_MAP.resize(m_FULLCALO_0p1x0p1_NETA, std::vector<double>(m_FULLCALO_0p1x0p1_NPHI, 0));
0350     m_FULLCALO_0p2x0p2_MAP.resize(m_FULLCALO_0p2x0p2_NETA, std::vector<double>(m_FULLCALO_0p2x0p2_NPHI, 0));
0351     m_FULLCALO_0p4x0p4_MAP.resize(m_FULLCALO_0p4x0p4_NETA, std::vector<double>(m_FULLCALO_0p4x0p4_NPHI, 0));
0352     m_FULLCALO_0p6x0p6_MAP.resize(m_FULLCALO_0p6x0p6_NETA, std::vector<double>(m_FULLCALO_0p6x0p6_NPHI, 0));
0353     m_FULLCALO_0p8x0p8_MAP.resize(m_FULLCALO_0p8x0p8_NETA, std::vector<double>(m_FULLCALO_0p8x0p8_NPHI, 0));
0354     m_FULLCALO_1p0x1p0_MAP.resize(m_FULLCALO_1p0x1p0_NETA, std::vector<double>(m_FULLCALO_1p0x1p0_NPHI, 0));
0355
0356     if (Verbosity() > 0)
0357     {
0358       std::cout << "CaloTriggerSim::process_event: determining phi range for 0.1x0.1 full calo map: " << m_FULLCALO_PHI_START << " to " << m_FULLCALO_PHI_END << std::endl;
0359       std::cout << "CaloTriggerSim::process_event: setting number of full calo window in eta / phi:" << std::endl;
0360       std::cout << "  0.1x0.1 are " << m_FULLCALO_0p1x0p1_NETA << " / " << m_FULLCALO_0p1x0p1_NPHI << ", ";
0361       std::cout << "0.2x0.2 are " << m_FULLCALO_0p2x0p2_NETA << " / " << m_FULLCALO_0p2x0p2_NPHI << ", ";
0362       std::cout << "0.4x0.4 are " << m_FULLCALO_0p4x0p4_NETA << " / " << m_FULLCALO_0p4x0p4_NPHI << ", ";
0363       std::cout << "0.6x0.6 are " << m_FULLCALO_0p6x0p6_NETA << " / " << m_FULLCALO_0p6x0p6_NPHI << ", ";
0364       std::cout << "0.8x0.8 are " << m_FULLCALO_0p8x0p8_NETA << " / " << m_FULLCALO_0p8x0p8_NPHI << ", ";
0365       std::cout << "1.0x1.0 are " << m_FULLCALO_1p0x1p0_NETA << " / " << m_FULLCALO_1p0x1p0_NPHI << std::endl;
0366     }
0367   }
0368
0369   // reset 0.1x0.1 map
0370   fill(m_FULLCALO_0p1x0p1_MAP.begin(), m_FULLCALO_0p1x0p1_MAP.end(), std::vector<double>(m_FULLCALO_0p1x0p1_NPHI, 0));
0371
0372   // iterate over EMCal towers, filling in the 0.1x0.1 region they contribute to
0373   for (unsigned int channel = 0; channel < ntowers_EM; channel++)
0374   {
0375     TowerInfo *tower = towersEM3->get_tower_at_channel(channel);
0376     unsigned int towerkey = towersEM3->encode_key(channel);
0377     int ieta = towersEM3->getTowerEtaBin(towerkey);
0378     int iphi = towersEM3->getTowerPhiBin(towerkey);
0379     const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0380     float this_phi = geomEM->get_tower_geometry(key)->get_phi();
0381     float this_eta = geomEM->get_tower_geometry(key)->get_eta();
0382
0383     if (this_phi < m_FULLCALO_PHI_START)
0384     {
0385       this_phi += 2 * M_PI;
0386     }
0387     if (this_phi > m_FULLCALO_PHI_END)
0388     {
0389       this_phi -= 2 * M_PI;
0390     }
0391
0392     // note: look up eta/phi index based on OHCal geometry, since this
0393     // defines the 0.1x0.1 regions
0394     int this_etabin = geomOH->get_etabin(this_eta);
0395     int this_phibin = geomOH->get_phibin(this_phi);
0396     double this_E = tower->get_energy();
0397
0398     m_FULLCALO_0p1x0p1_MAP[this_etabin][this_phibin] += this_E;
0399
0400     if (Verbosity() > 1 && tower->get_energy() > 1)
0401     {
0402       std::cout << "CaloTriggerSim::process_event: EMCal tower at eta / phi (added to fullcalo map with etabin / phibin ) / E = " << std::setprecision(6) << this_eta << " / " << this_phi << " ( " << this_etabin << " / " << this_phibin << " ) / " << this_E << std::endl;
0403     }
0404   }
0405
0406   // iterate over IHCal towers, filling in the 0.1x0.1 region they contribute to
0407
0408   unsigned int ntowers_IH = towersIH3->size();
0409   for (unsigned int channel = 0; channel < ntowers_IH; channel++)
0410   {
0411     TowerInfo *tower = towersIH3->get_tower_at_channel(channel);
0412     unsigned int towerkey = towersIH3->encode_key(channel);
0413     int ieta = towersIH3->getTowerEtaBin(towerkey);
0414     int iphi = towersIH3->getTowerPhiBin(towerkey);
0415     const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0416     float this_phi = geomIH->get_tower_geometry(key)->get_phi();
0417     float this_eta = geomIH->get_tower_geometry(key)->get_eta();
0418
0419     if (this_phi < m_FULLCALO_PHI_START)
0420     {
0421       this_phi += 2 * M_PI;
0422     }
0423     if (this_phi > m_FULLCALO_PHI_END)
0424     {
0425       this_phi -= 2 * M_PI;
0426     }
0427
0428     // note: look up eta/phi index based on OHCal geometry, even though I
0429     // think it is by construction the same as the IHCal geometry...
0430     int this_etabin = geomOH->get_etabin(this_eta);
0431     int this_phibin = geomOH->get_phibin(this_phi);
0432     double this_E = tower->get_energy();
0433
0434     m_FULLCALO_0p1x0p1_MAP[this_etabin][this_phibin] += this_E;
0435
0436     if (Verbosity() > 1 && tower->get_energy() > 0.5)
0437     {
0438       std::cout << "CaloTriggerSim::process_event: IHCal tower at eta / phi (added to fullcalo map with etabin / phibin ) / E = " << std::setprecision(6) << this_eta << " / " << this_phi << " ( " << this_etabin << " / " << this_phibin << " ) / " << this_E << std::endl;
0439     }
0440   }
0441
0442   // iterate over OHCal towers, filling in the 0.1x0.1 region they contribute to
0443
0444   unsigned int ntowers_OH = towersOH3->size();
0445   for (unsigned int channel = 0; channel < ntowers_OH; channel++)
0446   {
0447     TowerInfo *tower = towersOH3->get_tower_at_channel(channel);
0448     unsigned int towerkey = towersOH3->encode_key(channel);
0449     int ieta = towersOH3->getTowerEtaBin(towerkey);
0450     int iphi = towersOH3->getTowerPhiBin(towerkey);
0451     const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0452     float this_phi = geomIH->get_tower_geometry(key)->get_phi();
0453     float this_eta = geomIH->get_tower_geometry(key)->get_eta();
0454
0455     if (this_phi < m_FULLCALO_PHI_START)
0456     {
0457       this_phi += 2 * M_PI;
0458     }
0459     if (this_phi > m_FULLCALO_PHI_END)
0460     {
0461       this_phi -= 2 * M_PI;
0462     }
0463
0464     // note: use the nominal eta/phi index, since the fullcalo 0.1x0.1
0465     // map is defined by the OHCal geometry itself
0466     int this_etabin = geomOH->get_etabin(this_eta);
0467     int this_phibin = geomOH->get_phibin(this_phi);
0468     double this_E = tower->get_energy();
0469
0470     m_FULLCALO_0p1x0p1_MAP[this_etabin][this_phibin] += this_E;
0471
0472     if (Verbosity() > 1 && tower->get_energy() > 0.5)
0473     {
0474       std::cout << "CaloTriggerSim::process_event: OHCal tower at eta / phi (added to fullcalo map with etabin / phibin ) / E = " << std::setprecision(6) << this_eta << " / " << this_phi << " ( " << this_etabin << " / " << this_phibin << " ) / " << this_E << std::endl;
0475     }
0476   }
0477
0478   // reset 0.2x0.2 map and best
0479   fill(m_FULLCALO_0p2x0p2_MAP.begin(), m_FULLCALO_0p2x0p2_MAP.end(), std::vector<double>(m_FULLCALO_0p2x0p2_NPHI, 0));
0480
0481   m_FULLCALO_0p2x0p2_BEST_E = 0;
0482   m_FULLCALO_0p2x0p2_BEST_PHI = 0;
0483   m_FULLCALO_0p2x0p2_BEST_ETA = 0;
0484
0485   // now reconstruct (non-sliding) 0.2x0.2 map from 0.1x0.1 map
0486   for (std::vector<double>::size_type ieta = 0; ieta < m_FULLCALO_0p2x0p2_NETA; ieta++)
0487   {
0488     for (std::vector<double>::size_type iphi = 0; iphi < m_FULLCALO_0p2x0p2_NPHI; iphi++)
0489     {
0490       double this_sum = 0;
0491
0492       this_sum += m_FULLCALO_0p1x0p1_MAP[2 * ieta][2 * iphi];
0493       this_sum += m_FULLCALO_0p1x0p1_MAP[2 * ieta][(2 * iphi) + 1];  // 2 * iphi + 1 is safe, since m_FULLCALO_0p2x0p2_NPHI = m_FULLCALO_0p1x0p1_NPHI / 2
0494       this_sum += m_FULLCALO_0p1x0p1_MAP[(2 * ieta) + 1][2 * iphi];  // 2 * ieta + 1 is safe, since m_FULLCALO_0p2x0p2_NETA = m_FULLCALO_0p1x0p1_NETA / 2
0495       this_sum += m_FULLCALO_0p1x0p1_MAP[(2 * ieta) + 1][(2 * iphi) + 1];
0496
0497       // populate 0.2x0.2 map
0498       m_FULLCALO_0p2x0p2_MAP[ieta][iphi] = this_sum;
0499
0500       // to calculate the eta, phi position, take the average of that
0501       // of the contributing 0.1x0.1's (which are defined by the OHCal geometry)
0502       double this_eta = 0.5 * (geomOH->get_etacenter(2 * ieta) + geomOH->get_etacenter((2 * ieta) + 1));
0503       double this_phi = 0.5 * (geomOH->get_phicenter(2 * iphi) + geomOH->get_phicenter((2 * iphi) + 1));
0504
0505       if (Verbosity() > 1 && this_sum > 1)
0506       {
0507         std::cout << "CaloTriggerSim::process_event: FullCalo 0.2x0.2 window eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << ieta << " ) / " << this_phi << " ( " << iphi << " ) / " << this_sum << std::endl;
0508       }
0509
0510       if (this_sum > m_FULLCALO_0p2x0p2_BEST_E)
0511       {
0512         m_FULLCALO_0p2x0p2_BEST_E = this_sum;
0513         m_FULLCALO_0p2x0p2_BEST_PHI = this_phi;
0514         m_FULLCALO_0p2x0p2_BEST_ETA = this_eta;
0515       }
0516     }
0517   }
0518
0519   if (Verbosity() > 0)
0520   {
0521     std::cout << "CaloTriggerSim::process_event: best FullCalo 0.2x0.2 window is at eta / phi = " << m_FULLCALO_0p2x0p2_BEST_ETA << " / " << m_FULLCALO_0p2x0p2_BEST_PHI << " and E = " << m_FULLCALO_0p2x0p2_BEST_E << std::endl;
0522   }
0523
0524   // reset fullcalo 0.4x0.4 map & best
0525   fill(m_FULLCALO_0p4x0p4_MAP.begin(), m_FULLCALO_0p4x0p4_MAP.end(), std::vector<double>(m_FULLCALO_0p4x0p4_NPHI, 0));
0526
0527   m_FULLCALO_0p4x0p4_BEST_E = 0;
0528   m_FULLCALO_0p4x0p4_BEST_PHI = 0;
0529   m_FULLCALO_0p4x0p4_BEST_ETA = 0;
0530
0531   // now reconstruct (sliding) 0.4x0.4 map from 0.2x0.2 map
0532   for (int ieta = 0; ieta < m_FULLCALO_0p4x0p4_NETA; ieta++)
0533   {
0534     for (int iphi = 0; iphi < m_FULLCALO_0p4x0p4_NPHI; iphi++)
0535     {
0536       // for eta calculation, use position of corner tower and add 1.5
0537       // tower widths
0538       double this_eta = geomOH->get_etacenter(2 * ieta) + (1.5 * (geomOH->get_etacenter(1) - geomOH->get_etacenter(0)));
0539       // for phi calculation, use position of corner tower and add 1.5
0540       // tower widths
0541       double this_phi = geomOH->get_phicenter(2 * iphi) + (1.5 * (geomOH->get_phicenter(1) - geomOH->get_phicenter(0)));
0542
0543       double this_sum = 0;
0544
0545       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][iphi];
0546       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][iphi];  // 2 * ieta + 1 is safe, since m_FULLCALO_0p4x0p4_NETA = m_FULLCALO_0p4x0p4_NETA - 1
0547
0548       // add 1 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0549       // in case we have wrapped back around
0550       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0551       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0552
0553       m_FULLCALO_0p4x0p4_MAP[ieta][iphi] = this_sum;
0554
0555       if (Verbosity() > 1 && this_sum > 2)
0556       {
0557         std::cout << "CaloTriggerSim::process_event: FullCalo  0.4x0.4 tower eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << ieta << " ) / " << this_phi << " ( " << iphi << " ) / " << this_sum << std::endl;
0558       }
0559
0560       if (this_sum > m_FULLCALO_0p4x0p4_BEST_E)
0561       {
0562         m_FULLCALO_0p4x0p4_BEST_E = this_sum;
0563         m_FULLCALO_0p4x0p4_BEST_PHI = this_phi;
0564         m_FULLCALO_0p4x0p4_BEST_ETA = this_eta;
0565       }
0566     }
0567   }
0568
0569   if (Verbosity() > 0)
0570   {
0571     std::cout << "CaloTriggerSim::process_event: best FullCalo 0.4x0.4 window is at eta / phi = " << m_FULLCALO_0p4x0p4_BEST_ETA << " / " << m_FULLCALO_0p4x0p4_BEST_PHI << " and E = " << m_FULLCALO_0p4x0p4_BEST_E << std::endl;
0572   }
0573
0574   // reset fullcalo 0.6x0.6 map & best
0575   fill(m_FULLCALO_0p6x0p6_MAP.begin(), m_FULLCALO_0p6x0p6_MAP.end(), std::vector<double>(m_FULLCALO_0p6x0p6_NPHI, 0));
0576
0577   m_FULLCALO_0p6x0p6_BEST_E = 0;
0578   m_FULLCALO_0p6x0p6_BEST_PHI = 0;
0579   m_FULLCALO_0p6x0p6_BEST_ETA = 0;
0580
0581   // now reconstruct (sliding) 0.6x0.6 map from 0.2x0.2 map
0582   for (int ieta = 0; ieta < m_FULLCALO_0p6x0p6_NETA; ieta++)
0583   {
0584     for (int iphi = 0; iphi < m_FULLCALO_0p6x0p6_NPHI; iphi++)
0585     {
0586       // for eta calculation, use position of corner tower and add 2.5
0587       // tower widths
0588       double this_eta = geomOH->get_etacenter(2 * ieta) + (2.5 * (geomOH->get_etacenter(1) - geomOH->get_etacenter(0)));
0589       // for phi calculation, use position of corner tower and add 2.5
0590       // tower widths
0591       double this_phi = geomOH->get_phicenter(2 * iphi) + (2.5 * (geomOH->get_phicenter(1) - geomOH->get_phicenter(0)));
0592
0593       double this_sum = 0;
0594
0595       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][iphi];
0596       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][iphi];  // ieta + 1 is safe, since m_FULLCALO_0p6x0p6_NETA = m_FULLCALO_0p2x0p2_NETA - 2
0597       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][iphi];  // ieta + 2 is safe, since m_FULLCALO_0p6x0p6_NETA = m_FULLCALO_0p2x0p2_NETA - 2
0598
0599       // add 1 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0600       // in case we have wrapped back around
0601       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0602       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0603       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0604       // add 2 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0605       // in case we have wrapped back around
0606       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0607       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0608       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0609
0610       m_FULLCALO_0p6x0p6_MAP[ieta][iphi] = this_sum;
0611
0612       if (Verbosity() > 1 && this_sum > 3)
0613       {
0614         std::cout << "CaloTriggerSim::process_event: FullCalo  0.6x0.6 tower eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << ieta << " ) / " << this_phi << " ( " << iphi << " ) / " << this_sum << std::endl;
0615       }
0616
0617       if (this_sum > m_FULLCALO_0p6x0p6_BEST_E)
0618       {
0619         m_FULLCALO_0p6x0p6_BEST_E = this_sum;
0620         m_FULLCALO_0p6x0p6_BEST_PHI = this_phi;
0621         m_FULLCALO_0p6x0p6_BEST_ETA = this_eta;
0622       }
0623     }
0624   }
0625
0626   if (Verbosity() > 0)
0627   {
0628     std::cout << "CaloTriggerSim::process_event: best FullCalo 0.6x0.6 window is at eta / phi = " << m_FULLCALO_0p6x0p6_BEST_ETA << " / " << m_FULLCALO_0p6x0p6_BEST_PHI << " and E = " << m_FULLCALO_0p6x0p6_BEST_E << std::endl;
0629   }
0630
0631   // reset fullcalo 0.8x0.8 map & best
0632   fill(m_FULLCALO_0p8x0p8_MAP.begin(), m_FULLCALO_0p8x0p8_MAP.end(), std::vector<double>(m_FULLCALO_0p8x0p8_NPHI, 0));
0633
0634   m_FULLCALO_0p8x0p8_BEST_E = 0;
0635   m_FULLCALO_0p8x0p8_BEST_PHI = 0;
0636   m_FULLCALO_0p8x0p8_BEST_ETA = 0;
0637
0638   // now reconstruct (sliding) 0.8x0.8 map from 0.2x0.2 map
0639   for (int ieta = 0; ieta < m_FULLCALO_0p8x0p8_NETA; ieta++)
0640   {
0641     for (int iphi = 0; iphi < m_FULLCALO_0p8x0p8_NPHI; iphi++)
0642     {
0643       // for eta calculation, use position of corner tower and add 3.5
0644       // tower widths
0645       double this_eta = geomOH->get_etacenter(2 * ieta) + (3.5 * (geomOH->get_etacenter(1) - geomOH->get_etacenter(0)));
0646       // for phi calculation, use position of corner tower and add 3.5
0647       // tower widths
0648       double this_phi = geomOH->get_phicenter(2 * iphi) + (3.5 * (geomOH->get_phicenter(1) - geomOH->get_phicenter(0)));
0649
0650       double this_sum = 0;
0651
0652       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][iphi];
0653       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][iphi];  // ieta + 1 is safe, since m_FULLCALO_0p8x0p8_NETA = m_FULLCALO_0p2x0p2_NETA - 3
0654       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][iphi];  // ieta + 2 is safe, since m_FULLCALO_0p8x0p8_NETA = m_FULLCALO_0p2x0p2_NETA - 3
0655       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][iphi];  // ieta + 3 is safe, since m_FULLCALO_0p8x0p8_NETA = m_FULLCALO_0p2x0p2_NETA - 3
0656
0657       // add 1 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0658       // in case we have wrapped back around
0659       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0660       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0661       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0662       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0663       // add 2 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0664       // in case we have wrapped back around
0665       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0666       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0667       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0668       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0669       // add 3 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0670       // in case we have wrapped back around
0671       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0672       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0673       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0674       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0675
0676       m_FULLCALO_0p8x0p8_MAP[ieta][iphi] = this_sum;
0677
0678       if (Verbosity() > 1 && this_sum > 4)
0679       {
0680         std::cout << "CaloTriggerSim::process_event: FullCalo  0.8x0.8 tower eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << ieta << " ) / " << this_phi << " ( " << iphi << " ) / " << this_sum << std::endl;
0681       }
0682
0683       if (this_sum > m_FULLCALO_0p8x0p8_BEST_E)
0684       {
0685         m_FULLCALO_0p8x0p8_BEST_E = this_sum;
0686         m_FULLCALO_0p8x0p8_BEST_PHI = this_phi;
0687         m_FULLCALO_0p8x0p8_BEST_ETA = this_eta;
0688       }
0689     }
0690   }
0691
0692   if (Verbosity() > 0)
0693   {
0694     std::cout << "CaloTriggerSim::process_event: best FullCalo 0.8x0.8 window is at eta / phi = " << m_FULLCALO_0p8x0p8_BEST_ETA << " / " << m_FULLCALO_0p8x0p8_BEST_PHI << " and E = " << m_FULLCALO_0p8x0p8_BEST_E << std::endl;
0695   }
0696
0697   // reset fullcalo 1.0x1.0 map & best
0698   fill(m_FULLCALO_1p0x1p0_MAP.begin(), m_FULLCALO_1p0x1p0_MAP.end(), std::vector<double>(m_FULLCALO_1p0x1p0_NPHI, 0));
0699
0700   m_FULLCALO_1p0x1p0_BEST_E = 0;
0701   m_FULLCALO_1p0x1p0_BEST_PHI = 0;
0702   m_FULLCALO_1p0x1p0_BEST_ETA = 0;
0703
0704   // now reconstruct (sliding) 1.0x1.0 map from 0.2x0.2 map
0705   for (int ieta = 0; ieta < m_FULLCALO_1p0x1p0_NETA; ieta++)
0706   {
0707     for (int iphi = 0; iphi < m_FULLCALO_1p0x1p0_NPHI; iphi++)
0708     {
0709       // for eta calculation, use position of corner tower and add 4.5
0710       // tower widths
0711       double this_eta = geomOH->get_etacenter(2 * ieta) + (4.5 * (geomOH->get_etacenter(1) - geomOH->get_etacenter(0)));
0712       // for phi calculation, use position of corner tower and add 4.5
0713       // tower widths
0714       double this_phi = geomOH->get_phicenter(2 * iphi) + (4.5 * (geomOH->get_phicenter(1) - geomOH->get_phicenter(0)));
0715
0716       double this_sum = 0;
0717
0718       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][iphi];
0719       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][iphi];  // ieta + 1 is safe, since m_FULLCALO_1p0x1p0_NETA = m_FULLCALO_0p2x0p2_NETA - 4
0720       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][iphi];  // ieta + 2 is safe, since m_FULLCALO_1p0x1p0_NETA = m_FULLCALO_0p2x0p2_NETA - 4
0721       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][iphi];  // ieta + 3 is safe, since m_FULLCALO_1p0x1p0_NETA = m_FULLCALO_0p2x0p2_NETA - 4
0722       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 4][iphi];  // ieta + 4 is safe, since m_FULLCALO_1p0x1p0_NETA = m_FULLCALO_0p2x0p2_NETA - 4
0723
0724       // add 1 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0725       // in case we have wrapped back around
0726       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0727       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0728       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0729       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0730       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 4][(iphi + 1) % m_FULLCALO_0p2x0p2_NPHI];
0731       // add 2 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0732       // in case we have wrapped back around
0733       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0734       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0735       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0736       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0737       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 4][(iphi + 2) % m_FULLCALO_0p2x0p2_NPHI];
0738       // add 3 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0739       // in case we have wrapped back around
0740       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0741       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0742       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0743       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0744       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 4][(iphi + 3) % m_FULLCALO_0p2x0p2_NPHI];
0745       // add 4 to phi, but take modulus w.r.t. m_FULLCALO_0p2x0p2_NPHI
0746       // in case we have wrapped back around
0747       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta][(iphi + 4) % m_FULLCALO_0p2x0p2_NPHI];
0748       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 1][(iphi + 4) % m_FULLCALO_0p2x0p2_NPHI];
0749       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 2][(iphi + 4) % m_FULLCALO_0p2x0p2_NPHI];
0750       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 3][(iphi + 4) % m_FULLCALO_0p2x0p2_NPHI];
0751       this_sum += m_FULLCALO_0p2x0p2_MAP[ieta + 4][(iphi + 4) % m_FULLCALO_0p2x0p2_NPHI];
0752
0753       m_FULLCALO_1p0x1p0_MAP[ieta][iphi] = this_sum;
0754
0755       if (Verbosity() > 1 && this_sum > 5)
0756       {
0757         std::cout << "CaloTriggerSim::process_event: FullCalo  1.0x1.0 tower eta ( bin ) / phi ( bin ) / E = " << std::setprecision(6) << this_eta << " ( " << ieta << " ) / " << this_phi << " ( " << iphi << " ) / " << this_sum << std::endl;
0758       }
0759
0760       if (this_sum > m_FULLCALO_1p0x1p0_BEST_E)
0761       {
0762         m_FULLCALO_1p0x1p0_BEST_E = this_sum;
0763         m_FULLCALO_1p0x1p0_BEST_PHI = this_phi;
0764         m_FULLCALO_1p0x1p0_BEST_ETA = this_eta;
0765       }
0766     }
0767   }
0768
0769   if (Verbosity() > 0)
0770   {
0771     std::cout << "CaloTriggerSim::process_event: best FullCalo 1.0x1.0 window is at eta / phi = " << m_FULLCALO_1p0x1p0_BEST_ETA << " / " << m_FULLCALO_1p0x1p0_BEST_PHI << " and E = " << m_FULLCALO_1p0x1p0_BEST_E << std::endl;
0772   }
0773
0774   FillNode(topNode);
0775
0776   if (Verbosity() > 0)
0777   {
0778     std::cout << "CaloTriggerSim::process_event: exiting" << std::endl;
0779   }
0780
0781   return Fun4AllReturnCodes::EVENT_OK;
0782 }
0783
0784 int CaloTriggerSim::CreateNode(PHCompositeNode *topNode) const
0785 {
0786   PHNodeIterator iter(topNode);
0787
0788   // Looking for the DST node
0789   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0790   if (!dstNode)
0791   {
0792     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0793     return Fun4AllReturnCodes::ABORTRUN;
0794   }
0795
0796   // store the trigger stuff under a sub-node directory
0797   PHCompositeNode *trigNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "TRIG"));
0798   if (!trigNode)
0799   {
0800     trigNode = new PHCompositeNode("TRIG");
0801     dstNode->addNode(trigNode);
0802   }
0803
0804   // create the CaloTriggerInfo
0805   CaloTriggerInfo *triggerinfo = findNode::getClass<CaloTriggerInfo>(topNode, !m_EmulateTruncationFlag ? "CaloTriggerInfo" : "CaloTriggerInfo_Truncate");
0806   if (!triggerinfo)
0807   {
0808     triggerinfo = new CaloTriggerInfov1();
0809     PHIODataNode<PHObject> *TriggerNode = new PHIODataNode<PHObject>(triggerinfo, !m_EmulateTruncationFlag ? "CaloTriggerInfo" : "CaloTriggerInfo_Truncate", "PHObject");
0810     trigNode->addNode(TriggerNode);
0811   }
0812   else
0813   {
0814     std::cout << PHWHERE << "::ERROR - CaloTriggerInfo pre-exists, but should not" << std::endl;
0815     exit(-1);
0816   }
0817
0818   return Fun4AllReturnCodes::EVENT_OK;
0819 }
0820
0821 void CaloTriggerSim::FillNode(PHCompositeNode *topNode) const
0822 {
0823   CaloTriggerInfo *triggerinfo = findNode::getClass<CaloTriggerInfo>(topNode, !m_EmulateTruncationFlag ? "CaloTriggerInfo" : "CaloTriggerInfo_Truncate");
0824   if (!triggerinfo)
0825   {
0826     std::cout << " ERROR -- can't find CaloTriggerInfo node after it should have been created" << std::endl;
0827     return;
0828   }
0829
0830   triggerinfo->set_best_EMCal_2x2_E(m_EMCAL_2x2_BEST_E);
0831   triggerinfo->set_best_EMCal_2x2_eta(m_EMCAL_2x2_BEST_ETA);
0832   triggerinfo->set_best_EMCal_2x2_phi(m_EMCAL_2x2_BEST_PHI);
0833
0834   triggerinfo->set_best_EMCal_4x4_E(m_EMCAL_4x4_BEST_E);
0835   triggerinfo->set_best_EMCal_4x4_eta(m_EMCAL_4x4_BEST_ETA);
0836   triggerinfo->set_best_EMCal_4x4_phi(m_EMCAL_4x4_BEST_PHI);
0837
0838   triggerinfo->set_best2_EMCal_4x4_E(m_EMCAL_4x4_BEST2_E);
0839   triggerinfo->set_best2_EMCal_4x4_eta(m_EMCAL_4x4_BEST2_ETA);
0840   triggerinfo->set_best2_EMCal_4x4_phi(m_EMCAL_4x4_BEST2_PHI);
0841
0842   triggerinfo->set_best_FullCalo_0p2x0p2_E(m_FULLCALO_0p2x0p2_BEST_E);
0843   triggerinfo->set_best_FullCalo_0p2x0p2_eta(m_FULLCALO_0p2x0p2_BEST_ETA);
0844   triggerinfo->set_best_FullCalo_0p2x0p2_phi(m_FULLCALO_0p2x0p2_BEST_PHI);
0845
0846   triggerinfo->set_best_FullCalo_0p4x0p4_E(m_FULLCALO_0p4x0p4_BEST_E);
0847   triggerinfo->set_best_FullCalo_0p4x0p4_eta(m_FULLCALO_0p4x0p4_BEST_ETA);
0848   triggerinfo->set_best_FullCalo_0p4x0p4_phi(m_FULLCALO_0p4x0p4_BEST_PHI);
0849
0850   triggerinfo->set_best_FullCalo_0p6x0p6_E(m_FULLCALO_0p6x0p6_BEST_E);
0851   triggerinfo->set_best_FullCalo_0p6x0p6_eta(m_FULLCALO_0p6x0p6_BEST_ETA);
0852   triggerinfo->set_best_FullCalo_0p6x0p6_phi(m_FULLCALO_0p6x0p6_BEST_PHI);
0853
0854   triggerinfo->set_best_FullCalo_0p8x0p8_E(m_FULLCALO_0p8x0p8_BEST_E);
0855   triggerinfo->set_best_FullCalo_0p8x0p8_eta(m_FULLCALO_0p8x0p8_BEST_ETA);
0856   triggerinfo->set_best_FullCalo_0p8x0p8_phi(m_FULLCALO_0p8x0p8_BEST_PHI);
0857
0858   triggerinfo->set_best_FullCalo_1p0x1p0_E(m_FULLCALO_1p0x1p0_BEST_E);
0859   triggerinfo->set_best_FullCalo_1p0x1p0_eta(m_FULLCALO_1p0x1p0_BEST_ETA);
0860   triggerinfo->set_best_FullCalo_1p0x1p0_phi(m_FULLCALO_1p0x1p0_BEST_PHI);
0861
0862   return;
0863 }