3 #include <calobase/RawCluster.h>
4 #include <calobase/RawClusterContainer.h>
5 #include <calobase/RawTower.h>
6 #include <calobase/RawTowerContainer.h>
7 #include <calobase/RawTowerGeomContainer.h>
9 #include <phparameter/PHParameters.h>
35 :
SubsysReco(string(
"RawClusterPositionCorrection_") + name)
36 , _eclus_calib_params(string(
"eclus_params_") + name)
37 , _ecore_calib_params(string(
"ecore_params_") + name)
51 std::cout <<
"RawClusterPositionCorrection is running for clusters in the EMCal with eclus parameters:" << endl;
54 std::cout <<
"RawClusterPositionCorrection is running for clusters in the EMCal with ecore parameters:" << endl;
58 ostringstream paramname;
60 paramname <<
"number_of_bins";
68 for (
int j = 0; j <
bins; j++)
70 binvals.push_back(0. + j * 2. / (
float) (bins - 1));
73 for (
int i = 0; i < bins - 1; i++)
75 std::vector<double> dumvec;
77 for (
int j = 0; j < bins - 1; j++)
79 std::ostringstream calib_const_name;
80 calib_const_name.str(
"");
81 calib_const_name <<
"recalib_const_eta"
88 for (
int i = 0; i < bins - 1; i++)
90 std::vector<double> dumvec;
92 for (
int j = 0; j < bins - 1; j++)
94 std::ostringstream calib_const_name;
95 calib_const_name.str(
"");
96 calib_const_name <<
"recalib_const_eta"
110 std::cout <<
"Processing a NEW EVENT" << std::endl;
116 std::cout <<
"No " <<
_det_name <<
" Cluster Container found while in RawClusterPositionCorrection, can't proceed!!!" << std::endl;
123 std::cout <<
"No calibrated " <<
_det_name <<
" tower info found while in RawClusterPositionCorrection, can't proceed!!!" << std::endl;
127 string towergeomnodename =
"TOWERGEOM_" +
_det_name;
128 RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
131 cout <<
PHWHERE <<
": Could not find node " << towergeomnodename << endl;
140 for (iter = begin_end.first; iter != begin_end.second; ++iter)
147 std::vector<float> toweretas;
148 std::vector<float> towerphis;
149 std::vector<float> towerenergies;
155 for (toweriter = towers.first;
156 toweriter != towers.second;
166 toweretas.push_back(towereta);
167 towerphis.push_back(towerphi);
168 towerenergies.push_back(towerenergy);
171 int ntowers = toweretas.size();
172 assert(ntowers >= 1);
182 for (
int j = 0; j < ntowers; j++)
184 float energymult = towerenergies.at(j) * toweretas.at(j);
185 etamult += energymult;
186 etasum += towerenergies.at(j);
188 int phibin = towerphis.at(j);
190 if (phibin - towerphis.at(0) < -nphibin / 2.0)
192 else if (phibin - towerphis.at(0) > +nphibin / 2.0)
194 assert(
abs(phibin - towerphis.at(0)) <= nphibin / 2.0);
196 energymult = towerenergies.at(j) * phibin;
197 phimult += energymult;
198 phisum += towerenergies.at(j);
201 float avgphi = phimult / phisum;
202 float avgeta = etamult / etasum;
204 if (avgphi < 0) avgphi += nphibin;
207 float fmodphi = fmod(avgphi, 2.);
208 float fmodeta = fmod(avgeta, 2.);
215 for (
int j = 0; j <
bins - 1; j++)
219 for (
int j = 0; j < bins - 1; j++)
223 if ((phibin < 0 || etabin < 0) &&
Verbosity())
226 std::cout <<
"couldn't recalibrate cluster, something went wrong??" << std::endl;
229 float eclus_recalib_val = 1;
230 float ecore_recalib_val = 1;
231 if (phibin > -1 && etabin > -1)
237 assert(recalibcluster);
238 recalibcluster->
set_energy(clus_energy / eclus_recalib_val);
244 std::cout <<
"Input eclus cluster energy: " << clus_energy << endl;
245 std::cout <<
"Recalib value: " << eclus_recalib_val << endl;
246 std::cout <<
"Recalibrated eclus cluster energy: "
247 << clus_energy / eclus_recalib_val << endl;
248 std::cout <<
"Input ecore cluster energy: "
250 std::cout <<
"Recalib value: " << ecore_recalib_val << endl;
251 std::cout <<
"Recalibrated eclus cluster energy: "
252 << cluster->
get_ecore() / ecore_recalib_val << endl;
269 std::cerr <<
"DST Node missing, quitting" << std::endl;
270 throw std::runtime_error(
"failed to find DST node in RawClusterPositionCorrection::CreateNodeTree");
291 cemcNode->
addNode(clusterNode);
297 const string paramNodeName = string(
"eclus_Recalibration_" +
_det_name);
299 const string paramNodeName2 = string(
"ecore_Recalibration_" +
_det_name);
310 ostringstream param_name;
311 for (
int i = 0; i <
bins - 1; i++)
313 for (
int j = 0; j < bins - 1; j++)
316 param_name <<
"recalib_const_eta"