forked from sPHENIX-Collaboration/coresoftware
/
RawClusterPositionCorrection.cc
323 lines (266 loc) · 10.1 KB
/
RawClusterPositionCorrection.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
#include "RawClusterPositionCorrection.h"
#include <calobase/RawCluster.h>
#include <calobase/RawClusterContainer.h>
#include <calobase/RawTower.h>
#include <calobase/RawTowerContainer.h>
#include <calobase/RawTowerGeomContainer.h>
#include <phparameter/PHParameters.h>
#include <fun4all/Fun4AllReturnCodes.h>
#include <fun4all/SubsysReco.h>
#include <phool/getClass.h>
#include <phool/PHCompositeNode.h>
#include <phool/PHIODataNode.h>
#include <phool/PHNode.h>
#include <phool/PHNodeIterator.h>
#include <phool/PHObject.h>
#include <phool/phool.h>
#include <cassert>
#include <fstream>
#include <iostream>
#include <map>
#include <memory>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility> // for pair
using namespace std;
RawClusterPositionCorrection::RawClusterPositionCorrection(const std::string &name)
: SubsysReco(string("RawClusterPositionCorrection_") + name)
, _eclus_calib_params(string("eclus_params_") + name)
, _ecore_calib_params(string("ecore_params_") + name)
, _det_name(name)
, bins(17) //default bins to be 17 to set default recalib parameters to 1
{
SetDefaultParameters(_eclus_calib_params);
SetDefaultParameters(_ecore_calib_params);
}
int RawClusterPositionCorrection::InitRun(PHCompositeNode *topNode)
{
CreateNodeTree(topNode);
if (Verbosity())
{
std::cout << "RawClusterPositionCorrection is running for clusters in the EMCal with eclus parameters:" << endl;
_eclus_calib_params.Print();
std::cout << "RawClusterPositionCorrection is running for clusters in the EMCal with ecore parameters:" << endl;
_ecore_calib_params.Print();
}
//now get the actual number of bins in the calib file
ostringstream paramname;
paramname.str("");
paramname << "number_of_bins";
//+1 because I use bins as the total number of bin boundaries
//i.e. 16 bins corresponds to 17 bin boundaries
bins = _eclus_calib_params.get_int_param(paramname.str()) + 1;
//set bin boundaries
for (int j = 0; j < bins; j++)
{
binvals.push_back(0. + j * 2. / (float) (bins - 1));
}
for (int i = 0; i < bins - 1; i++)
{
std::vector<double> dumvec;
for (int j = 0; j < bins - 1; j++)
{
std::ostringstream calib_const_name;
calib_const_name.str("");
calib_const_name << "recalib_const_eta"
<< i << "_phi" << j;
dumvec.push_back(_eclus_calib_params.get_double_param(calib_const_name.str()));
}
eclus_calib_constants.push_back(dumvec);
}
for (int i = 0; i < bins - 1; i++)
{
std::vector<double> dumvec;
for (int j = 0; j < bins - 1; j++)
{
std::ostringstream calib_const_name;
calib_const_name.str("");
calib_const_name << "recalib_const_eta"
<< i << "_phi" << j;
dumvec.push_back(_ecore_calib_params.get_double_param(calib_const_name.str()));
}
ecore_calib_constants.push_back(dumvec);
}
return Fun4AllReturnCodes::EVENT_OK;
}
int RawClusterPositionCorrection::process_event(PHCompositeNode *topNode)
{
if (Verbosity())
{
std::cout << "Processing a NEW EVENT" << std::endl;
}
RawClusterContainer *rawclusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_" + _det_name);
if (!rawclusters)
{
std::cout << "No " << _det_name << " Cluster Container found while in RawClusterPositionCorrection, can't proceed!!!" << std::endl;
return Fun4AllReturnCodes::ABORTEVENT;
}
RawTowerContainer *_towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_" + _det_name);
if (!_towers)
{
std::cout << "No calibrated " << _det_name << " tower info found while in RawClusterPositionCorrection, can't proceed!!!" << std::endl;
return Fun4AllReturnCodes::ABORTEVENT;
}
string towergeomnodename = "TOWERGEOM_" + _det_name;
RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
if (!towergeom)
{
cout << PHWHERE << ": Could not find node " << towergeomnodename << endl;
return Fun4AllReturnCodes::ABORTEVENT;
}
const int nphibin = towergeom->get_phibins();
//loop over the clusters
RawClusterContainer::ConstRange begin_end = rawclusters->getClusters();
RawClusterContainer::ConstIterator iter;
for (iter = begin_end.first; iter != begin_end.second; ++iter)
{
// RawClusterDefs::keytype key = iter->first;
RawCluster *cluster = iter->second;
float clus_energy = cluster->get_energy();
std::vector<float> toweretas;
std::vector<float> towerphis;
std::vector<float> towerenergies;
//loop over the towers in the cluster
RawCluster::TowerConstRange towers = cluster->get_towers();
RawCluster::TowerConstIterator toweriter;
for (toweriter = towers.first;
toweriter != towers.second;
++toweriter)
{
RawTower *tower = _towers->getTower(toweriter->first);
int towereta = tower->get_bineta();
int towerphi = tower->get_binphi();
double towerenergy = tower->get_energy();
//put the etabin, phibin, and energy into the corresponding vectors
toweretas.push_back(towereta);
towerphis.push_back(towerphi);
towerenergies.push_back(towerenergy);
}
int ntowers = toweretas.size();
assert(ntowers >= 1);
//loop over the towers to determine the energy
//weighted eta and phi position of the cluster
float etamult = 0;
float etasum = 0;
float phimult = 0;
float phisum = 0;
for (int j = 0; j < ntowers; j++)
{
float energymult = towerenergies.at(j) * toweretas.at(j);
etamult += energymult;
etasum += towerenergies.at(j);
int phibin = towerphis.at(j);
if (phibin - towerphis.at(0) < -nphibin / 2.0)
phibin += nphibin;
else if (phibin - towerphis.at(0) > +nphibin / 2.0)
phibin -= nphibin;
assert(abs(phibin - towerphis.at(0)) <= nphibin / 2.0);
energymult = towerenergies.at(j) * phibin;
phimult += energymult;
phisum += towerenergies.at(j);
}
float avgphi = phimult / phisum;
float avgeta = etamult / etasum;
if (avgphi < 0) avgphi += nphibin;
//this determines the position of the cluster in the 2x2 block
float fmodphi = fmod(avgphi, 2.);
float fmodeta = fmod(avgeta, 2.);
//determine the bin number
//2 is here since we divide the 2x2 block into 16 bins in eta/phi
int etabin = -99;
int phibin = -99;
for (int j = 0; j < bins - 1; j++)
if (fmodphi >= binvals.at(j) && fmodphi <= binvals.at(j + 1))
phibin = j;
for (int j = 0; j < bins - 1; j++)
if (fmodeta >= binvals.at(j) && fmodeta <= binvals.at(j + 1))
etabin = j;
if ((phibin < 0 || etabin < 0) && Verbosity())
{
if (Verbosity())
std::cout << "couldn't recalibrate cluster, something went wrong??" << std::endl;
}
float eclus_recalib_val = 1;
float ecore_recalib_val = 1;
if (phibin > -1 && etabin > -1)
{
eclus_recalib_val = eclus_calib_constants.at(etabin).at(phibin);
ecore_recalib_val = ecore_calib_constants.at(etabin).at(phibin);
}
RawCluster *recalibcluster = dynamic_cast<RawCluster *>(cluster->CloneMe());
assert(recalibcluster);
recalibcluster->set_energy(clus_energy / eclus_recalib_val);
recalibcluster->set_ecore(cluster->get_ecore() / ecore_recalib_val);
_recalib_clusters->AddCluster(recalibcluster);
if (Verbosity() && clus_energy > 1)
{
std::cout << "Input eclus cluster energy: " << clus_energy << endl;
std::cout << "Recalib value: " << eclus_recalib_val << endl;
std::cout << "Recalibrated eclus cluster energy: "
<< clus_energy / eclus_recalib_val << endl;
std::cout << "Input ecore cluster energy: "
<< cluster->get_ecore() << endl;
std::cout << "Recalib value: " << ecore_recalib_val << endl;
std::cout << "Recalibrated eclus cluster energy: "
<< cluster->get_ecore() / ecore_recalib_val << endl;
}
}
return Fun4AllReturnCodes::EVENT_OK;
}
void RawClusterPositionCorrection::CreateNodeTree(PHCompositeNode *topNode)
{
PHNodeIterator iter(topNode);
//Get the DST Node
PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
//Check that it is there
if (!dstNode)
{
std::cerr << "DST Node missing, quitting" << std::endl;
throw std::runtime_error("failed to find DST node in RawClusterPositionCorrection::CreateNodeTree");
}
//Get the _det_name subnode
PHCompositeNode *cemcNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", _det_name));
//Check that it is there
if (!cemcNode)
{
cemcNode = new PHCompositeNode(_det_name);
dstNode->addNode(cemcNode);
}
//Check to see if the cluster recalib node is on the nodetree
_recalib_clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_RECALIB_" + _det_name);
//If not, make it and add it to the _det_name subnode
if (!_recalib_clusters)
{
_recalib_clusters = new RawClusterContainer();
PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_recalib_clusters, "CLUSTER_POS_COR_" + _det_name, "PHObject");
cemcNode->addNode(clusterNode);
}
//put the recalib parameters on the node tree
PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
assert(parNode);
const string paramNodeName = string("eclus_Recalibration_" + _det_name);
_eclus_calib_params.SaveToNodeTree(parNode, paramNodeName);
const string paramNodeName2 = string("ecore_Recalibration_" + _det_name);
_ecore_calib_params.SaveToNodeTree(parNode, paramNodeName2);
}
int RawClusterPositionCorrection::End(PHCompositeNode */*topNode*/)
{
return Fun4AllReturnCodes::EVENT_OK;
}
void RawClusterPositionCorrection::SetDefaultParameters(PHParameters ¶m)
{
param.set_int_param("number_of_bins", 17);
ostringstream param_name;
for (int i = 0; i < bins - 1; i++)
{
for (int j = 0; j < bins - 1; j++)
{
param_name.str("");
param_name << "recalib_const_eta"
<< i << "_phi" << j;
//default to 1, i.e. no recalibration
param.set_double_param(param_name.str(), 1.0);
}
}
}