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