28 :
SubsysReco(name), _centrality_calibration_params(name)
35 std::cout <<
"PHG4CentralityReco::InitRun : enter " << std::endl;
41 catch (std::exception &
e)
43 std::cout <<
PHWHERE <<
": " << e.what() << std::endl;
47 auto bhits = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_BBC");
49 std::cout <<
"PHG4CentralityReco::InitRun : cannot find G4HIT_BBC, will not use MBD centrality" << std::endl;
51 auto ehits = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_EPD");
53 std::cout <<
"PHG4CentralityReco::InitRun : cannot find G4HIT_EPD, will not use sEPD centrality" << std::endl;
57 std::cout <<
"PHG4CentralityReco::InitRun : Centrality calibration description : " << std::endl <<
" ";
61 for (
int n = 0;
n < 101;
n++) {
62 std::ostringstream
s1; s1 <<
"epd_centile_" <<
n;
69 for (
int n = 0;
n < 101;
n++) {
70 std::ostringstream s2; s2 <<
"mbd_centile_" <<
n;
77 for (
int n = 0;
n < 101;
n++) {
78 std::ostringstream s3; s3 <<
"bimp_centile_" <<
n;
88 std::cout <<
"PHG4CentralityReco::InitRun : no centrality calibration found!" << std::endl;
99 auto event_header = findNode::getClass<EventHeaderv1>(topNode,
"EventHeader" );
100 if ( event_header ) {
101 _bimp = event_header->get_floatval(
"bimp");
104 std::cout <<
"PHG4CentralityReco::process_event : Hijing impact parameter b = " <<
_bimp << std::endl;
109 std::cout <<
"PHG4CentralityReco::process_event : No Hijing impact parameter info, setting b = 101" << std::endl;
117 auto bhits = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_BBC");
121 auto brange = bhits->getHits();
122 for (
auto it = brange.first;
it != brange.second; ++
it)
124 if ((
it->second->get_t(0) > -50) && (
it->second->get_t(1) < 50))
127 int id =
it->second->get_layer();
128 if ((
id & 0x40) == 0)
136 std::cout <<
"PHG4CentralityReco::process_event : MBD Sum Charge N / S / N+S = " <<
_mbd_N <<
" / " <<
_mbd_S <<
" / " <<
_mbd_NS << std::endl;
141 std::cout <<
"PHG4CentralityReco::process_event : No MBD info, setting all Sum Charges = 0" << std::endl;
148 auto ehits = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_EPD");
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))
157 int id =
it->second->get_scint_id();
158 if ((
id & 0x200) == 0)
165 std::cout <<
"PHG4CentralityReco::process_event : sEPD Sum Energy N / S / N+S = " <<
_epd_N <<
" / " <<
_epd_S <<
" / " <<
_epd_NS << std::endl;
170 std::cout <<
"PHG4CentralityReco::process_event : No sEPD info, setting all Sum Energies = 0" << std::endl;
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;
179 float low_epd_val = -10000;
180 float high_epd_val = 10000;
181 int low_epd_centile = -1;
182 int high_epd_centile = -1;
185 float signal =
it->first;
186 int cent =
it->second;
188 if (signal < _epd_NS && signal > low_epd_val) {
189 low_epd_val = signal;
190 low_epd_centile = cent;
192 if (signal >
_epd_NS && signal < high_epd_val) {
193 high_epd_val = signal;
194 high_epd_centile = cent;
199 if ( low_epd_centile >= 0 && high_epd_centile >= 0 ) {
200 _epd_cent = ( low_epd_centile + high_epd_centile ) / 2.0;
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;
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;
210 float low_mbd_val = -10000;
211 float high_mbd_val = 10000;
212 int low_mbd_centile = -1;
213 int high_mbd_centile = -1;
216 float signal =
it->first;
217 int cent =
it->second;
219 if (signal < _mbd_NS && signal > low_mbd_val) {
220 low_mbd_val = signal;
221 low_mbd_centile = cent;
223 if (signal >
_mbd_NS && signal < high_mbd_val) {
224 high_mbd_val = signal;
225 high_mbd_centile = cent;
230 if ( low_mbd_centile >= 0 && high_mbd_centile >= 0 ) {
231 _mbd_cent = ( low_mbd_centile + high_mbd_centile ) / 2.0;
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;
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;
241 float low_bimp_val = -1;
242 float high_bimp_val = 10000;
243 int low_bimp_centile = -1;
244 int high_bimp_centile = -1;
247 float signal =
it->first;
248 int cent =
it->second;
250 if (signal < _bimp && signal > low_bimp_val) {
251 low_bimp_val = signal;
252 low_bimp_centile = cent;
254 if (signal >
_bimp && signal < high_bimp_val) {
255 high_bimp_val = signal;
256 high_bimp_centile = cent;
261 if ( low_bimp_centile >= 0 && high_bimp_centile >= 0 ) {
262 _bimp_cent = ( low_bimp_centile + high_bimp_centile ) / 2.0;
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;
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;
286 CentralityInfo *cent = findNode::getClass<CentralityInfo>(topNode,
"CentralityInfo");
289 std::cout <<
" ERROR -- can't find CentralityInfo node after it should have been created" << std::endl;
315 std::cerr <<
PHWHERE <<
"DST Node missing, doing nothing." << std::endl;
316 throw std::runtime_error(
"Failed to find DST node in PHG4CentralityReco::CreateNode");