ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4CentralityReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CentralityReco.cc
1 #include "PHG4CentralityReco.h"
2 
3 #include <centrality/CentralityInfo.h> // for CentralityInfo, CentralityI...
4 #include <centrality/CentralityInfov1.h> // for CentralityInfov1
5 
7 #include <fun4all/SubsysReco.h>
8 
10 #include <phool/PHIODataNode.h>
11 #include <phool/PHNode.h>
12 #include <phool/PHNodeIterator.h>
13 #include <phool/PHObject.h> // for PHObject
14 #include <phool/getClass.h>
15 #include <phool/phool.h> // for PHWHERE
16 
18 #include <g4main/PHG4Hit.h>
20 
21 #include <iostream> // for operator<<, basic_ostream
22 #include <map> // for _Rb_tree_const_iterator
23 #include <stdexcept> // for runtime_error
24 #include <utility> // for pair
25 #include <sstream>
26 
28  : SubsysReco(name), _centrality_calibration_params(name)
29 {
30 }
31 
33 {
34  if (Verbosity() >= 1)
35  std::cout << "PHG4CentralityReco::InitRun : enter " << std::endl;
36 
37  try
38  {
39  CreateNode(topNode);
40  }
41  catch (std::exception &e)
42  {
43  std::cout << PHWHERE << ": " << e.what() << std::endl;
44  throw;
45  }
46 
47  auto bhits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_BBC");
48  if (!bhits)
49  std::cout << "PHG4CentralityReco::InitRun : cannot find G4HIT_BBC, will not use MBD centrality" << std::endl;
50 
51  auto ehits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_EPD");
52  if (!ehits)
53  std::cout << "PHG4CentralityReco::InitRun : cannot find G4HIT_EPD, will not use sEPD centrality" << std::endl;
54 
56  if (Verbosity() >= 1 ) {
57  std::cout << "PHG4CentralityReco::InitRun : Centrality calibration description : " << std::endl << " ";
58  std::cout << _centrality_calibration_params.get_string_param("description") << std::endl;
59 
60  // search for possible centile definitions
61  for (int n = 0; n < 101; n++) {
62  std::ostringstream s1; s1 << "epd_centile_" << n;
63  if ( _centrality_calibration_params.exist_double_param( s1.str().c_str() ) ) {
65  if (Verbosity() >= 2 )
66  std::cout << "PHG4CentralityReco::InitRun : sEPD centrality calibration, centile " << n << "% is " << _centrality_calibration_params.get_double_param( s1.str().c_str() ) << std::endl;
67  }
68  }
69  for (int n = 0; n < 101; n++) {
70  std::ostringstream s2; s2 << "mbd_centile_" << n;
71  if ( _centrality_calibration_params.exist_double_param( s2.str().c_str() ) ) {
73  if (Verbosity() >= 2 )
74  std::cout << "PHG4CentralityReco::InitRun : MBD centrality calibration, centile " << n << "% is " << _centrality_calibration_params.get_double_param( s2.str().c_str() ) << std::endl;
75  }
76  }
77  for (int n = 0; n < 101; n++) {
78  std::ostringstream s3; s3 << "bimp_centile_" << n;
79  if ( _centrality_calibration_params.exist_double_param( s3.str().c_str() ) ) {
81  if (Verbosity() >= 2 )
82  std::cout << "PHG4CentralityReco::InitRun : b (impact parameter) centrality calibration, centile " << n << "% is " << _centrality_calibration_params.get_double_param( s3.str().c_str() ) << std::endl;
83  }
84  }
85 
86  }
87  } else {
88  std::cout << "PHG4CentralityReco::InitRun : no centrality calibration found!" << std::endl;
89  }
90 
92 }
93 
95 {
96 
97  _bimp = 101;
98 
99  auto event_header = findNode::getClass<EventHeaderv1>(topNode, "EventHeader" );
100  if ( event_header ) {
101  _bimp = event_header->get_floatval("bimp");
102 
103  if (Verbosity() >= 5)
104  std::cout << "PHG4CentralityReco::process_event : Hijing impact parameter b = " << _bimp << std::endl;
105  }
106  else
107  {
108  if (Verbosity() >= 5)
109  std::cout << "PHG4CentralityReco::process_event : No Hijing impact parameter info, setting b = 101" << std::endl;
110  }
111 
112 
113  _mbd_N = 0;
114  _mbd_S = 0;
115  _mbd_NS = 0;
116 
117  auto bhits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_BBC");
118 
119  if (bhits)
120  {
121  auto brange = bhits->getHits();
122  for (auto it = brange.first; it != brange.second; ++it)
123  {
124  if ((it->second->get_t(0) > -50) && (it->second->get_t(1) < 50))
125  {
126  _mbd_NS += it->second->get_edep();
127  int id = it->second->get_layer();
128  if ((id & 0x40) == 0)
129  _mbd_N += it->second->get_edep();
130  else
131  _mbd_S += it->second->get_edep();
132  }
133  }
134 
135  if (Verbosity() >= 5)
136  std::cout << "PHG4CentralityReco::process_event : MBD Sum Charge N / S / N+S = " << _mbd_N << " / " << _mbd_S << " / " << _mbd_NS << std::endl;
137  }
138  else
139  {
140  if (Verbosity() >= 5)
141  std::cout << "PHG4CentralityReco::process_event : No MBD info, setting all Sum Charges = 0" << std::endl;
142  }
143 
144  _epd_N = 0;
145  _epd_S = 0;
146  _epd_NS = 0;
147 
148  auto ehits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_EPD");
149 
150  if (ehits)
151  {
152  auto erange = ehits->getHits();
153  for (auto it = erange.first; it != erange.second; ++it)
154  if ((it->second->get_t(0) > -50) && (it->second->get_t(1) < 50))
155  {
156  _epd_NS += it->second->get_edep();
157  int id = it->second->get_scint_id();
158  if ((id & 0x200) == 0)
159  _epd_N += it->second->get_edep();
160  else
161  _epd_S += it->second->get_edep();
162  }
163 
164  if (Verbosity() >= 5)
165  std::cout << "PHG4CentralityReco::process_event : sEPD Sum Energy N / S / N+S = " << _epd_N << " / " << _epd_S << " / " << _epd_NS << std::endl;
166  }
167  else
168  {
169  if (Verbosity() >= 5)
170  std::cout << "PHG4CentralityReco::process_event : No sEPD info, setting all Sum Energies = 0" << std::endl;
171  }
172 
173  if (Verbosity() >= 1)
174  std::cout << "PHG4CentralityReco::process_event : summary MBD (N, S, N+S) = (" << _mbd_N << ", " << _mbd_S << ", " << _mbd_NS << "), sEPD (N, S, N+S) = (" << _epd_N << ", " << _epd_S << ", " << _epd_NS << ")" << std::endl;
175 
177 
178  // sEPD centrality
179  float low_epd_val = -10000;
180  float high_epd_val = 10000;
181  int low_epd_centile = -1;
182  int high_epd_centile = -1;
183 
184  for(std::map<float,int>::iterator it = _cent_cal_epd.begin(); it != _cent_cal_epd.end(); ++it) {
185  float signal = it->first;
186  int cent = it->second;
187 
188  if (signal < _epd_NS && signal > low_epd_val) {
189  low_epd_val = signal;
190  low_epd_centile = cent;
191  }
192  if (signal > _epd_NS && signal < high_epd_val) {
193  high_epd_val = signal;
194  high_epd_centile = cent;
195  }
196 
197  } // close iterate through sEPD cuts
198 
199  if ( low_epd_centile >= 0 && high_epd_centile >= 0 ) {
200  _epd_cent = ( low_epd_centile + high_epd_centile ) / 2.0;
201  if (Verbosity() >= 10)
202  std::cout << "PHG4CentralityReco::process_event : lower EPD value is " << low_epd_val << " (" << low_epd_centile << "%), higher is " << high_epd_val << " (" << high_epd_centile << "%), assigning " << _epd_cent << "%" << std::endl;
203  } else {
204  _epd_cent = 101;
205  if (Verbosity() >= 5)
206  std::cout << "PHG4CentralityReco::process_event : not able to map EPD value to a centrality. debug info = " << low_epd_val << "/" << low_epd_centile << "/" << high_epd_val << "/" << high_epd_centile << std::endl;
207  }
208 
209  // MBD centrality
210  float low_mbd_val = -10000;
211  float high_mbd_val = 10000;
212  int low_mbd_centile = -1;
213  int high_mbd_centile = -1;
214 
215  for(std::map<float,int>::iterator it = _cent_cal_mbd.begin(); it != _cent_cal_mbd.end(); ++it) {
216  float signal = it->first;
217  int cent = it->second;
218 
219  if (signal < _mbd_NS && signal > low_mbd_val) {
220  low_mbd_val = signal;
221  low_mbd_centile = cent;
222  }
223  if (signal > _mbd_NS && signal < high_mbd_val) {
224  high_mbd_val = signal;
225  high_mbd_centile = cent;
226  }
227 
228  } // close iterate through MBD cuts
229 
230  if ( low_mbd_centile >= 0 && high_mbd_centile >= 0 ) {
231  _mbd_cent = ( low_mbd_centile + high_mbd_centile ) / 2.0;
232  if (Verbosity() >= 10)
233  std::cout << "PHG4CentralityReco::process_event : lower MBD value is " << low_mbd_val << " (" << low_mbd_centile << "%), higher is " << high_mbd_val << " (" << high_mbd_centile << "%), assigning " << _mbd_cent << "%" << std::endl;
234  } else {
235  _mbd_cent = 101;
236  if (Verbosity() >= 5)
237  std::cout << "PHG4CentralityReco::process_event : not able to map MBD value to a centrality. debug info = " << low_mbd_val << "/" << low_mbd_centile << "/" << high_mbd_val << "/" << high_mbd_centile << std::endl;
238  }
239 
240  // b (impact parameter) centrality
241  float low_bimp_val = -1;
242  float high_bimp_val = 10000;
243  int low_bimp_centile = -1;
244  int high_bimp_centile = -1;
245 
246  for(std::map<float,int>::iterator it = _cent_cal_bimp.begin(); it != _cent_cal_bimp.end(); ++it) {
247  float signal = it->first;
248  int cent = it->second;
249 
250  if (signal < _bimp && signal > low_bimp_val) {
251  low_bimp_val = signal;
252  low_bimp_centile = cent;
253  }
254  if (signal > _bimp && signal < high_bimp_val) {
255  high_bimp_val = signal;
256  high_bimp_centile = cent;
257  }
258 
259  } // close iterate through bimp cuts
260 
261  if ( low_bimp_centile >= 0 && high_bimp_centile >= 0 ) {
262  _bimp_cent = ( low_bimp_centile + high_bimp_centile ) / 2.0;
263  if (Verbosity() >= 10)
264  std::cout << "PHG4CentralityReco::process_event : lower b value is " << low_bimp_val << " (" << low_bimp_centile << "%), higher is " << high_bimp_val << " (" << high_bimp_centile << "%), assigning " << _bimp_cent << "%" << std::endl;
265  } else {
266  _bimp_cent = 101;
267  if (Verbosity() >= 5)
268  std::cout << "PHG4CentralityReco::process_event : not able to map b value to a centrality. debug info = " << low_bimp_val << "/" << low_bimp_centile << "/" << high_bimp_val << "/" << high_bimp_centile << std::endl;
269  }
270 
271  } // close centrality calibration
272 
273 
274  FillNode(topNode);
275 
277 }
278 
280 {
282 }
283 
285 {
286  CentralityInfo *cent = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
287  if (!cent)
288  {
289  std::cout << " ERROR -- can't find CentralityInfo node after it should have been created" << std::endl;
290  return;
291  }
292  else
293  {
294  cent->set_quantity(CentralityInfo::PROP::mbd_N, _mbd_N);
295  cent->set_quantity(CentralityInfo::PROP::mbd_S, _mbd_S);
296  cent->set_quantity(CentralityInfo::PROP::mbd_NS, _mbd_NS);
297  cent->set_quantity(CentralityInfo::PROP::epd_N, _epd_N);
298  cent->set_quantity(CentralityInfo::PROP::epd_S, _epd_S);
299  cent->set_quantity(CentralityInfo::PROP::epd_NS, _epd_NS);
300  cent->set_quantity(CentralityInfo::PROP::bimp, _bimp);
301 
302  cent->set_centile(CentralityInfo::PROP::epd_NS, _epd_cent);
303  cent->set_centile(CentralityInfo::PROP::mbd_NS, _mbd_cent);
304  cent->set_centile(CentralityInfo::PROP::bimp, _bimp_cent);
305  }
306 }
307 
309 {
310  PHNodeIterator iter(topNode);
311 
312  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
313  if (!dstNode)
314  {
315  std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
316  throw std::runtime_error("Failed to find DST node in PHG4CentralityReco::CreateNode");
317  }
318 
319  PHNodeIterator dstiter(dstNode);
320  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "GLOBAL"));
321  if (!DetNode)
322  {
323  DetNode = new PHCompositeNode("GLOBAL");
324  dstNode->addNode(DetNode);
325  }
326 
327  CentralityInfo *cent = new CentralityInfov1();
328 
329  PHIODataNode<PHObject> *centNode = new PHIODataNode<PHObject>(cent, "CentralityInfo", "PHObject");
330  DetNode->addNode(centNode);
331 }