g4simulation/g4centrality/PHG4CentralityReco.cc

0001 #include "PHG4CentralityReco.h"
0002
0003 #include <centrality/CentralityInfo.h>    // for CentralityInfo, CentralityI...
0004 #include <centrality/CentralityInfov1.h>  // for CentralityInfov1
0005
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/SubsysReco.h>
0008
0009 #include <phool/PHCompositeNode.h>
0010 #include <phool/PHIODataNode.h>
0011 #include <phool/PHNode.h>
0012 #include <phool/PHNodeIterator.h>
0013 #include <phool/PHObject.h>  // for PHObject
0014 #include <phool/getClass.h>
0015 #include <phool/phool.h>  // for PHWHERE
0016
0017 #include <epd/EPDDefs.h>
0018 #include <ffaobjects/EventHeaderv1.h>
0019 #include <g4main/PHG4Hit.h>
0020 #include <g4main/PHG4HitContainer.h>
0021
0022 #include <iostream>  // for operator<<, basic_ostream
0023 #include <map>       // for _Rb_tree_const_iterator
0024 #include <sstream>
0025 #include <stdexcept>  // for runtime_error
0026 #include <utility>    // for pair
0027
0028 PHG4CentralityReco::PHG4CentralityReco(const std::string &name)
0029   : SubsysReco(name)
0030   , _centrality_calibration_params(name)
0031 {
0032 }
0033
0034 int PHG4CentralityReco::InitRun(PHCompositeNode *topNode)
0035 {
0036   if (Verbosity() >= 1)
0037   {
0038     std::cout << "PHG4CentralityReco::InitRun : enter " << std::endl;
0039   }
0040
0041   try
0042   {
0043     CreateNode(topNode);
0044   }
0045   catch (std::exception &e)
0046   {
0047     std::cout << PHWHERE << ": " << e.what() << std::endl;
0048     throw;
0049   }
0050
0051   auto bhits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_BBC");
0052   if (!bhits)
0053   {
0054     std::cout << "PHG4CentralityReco::InitRun : cannot find G4HIT_BBC, will not use MBD centrality" << std::endl;
0055   }
0056
0057   auto ehits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_EPD");
0058   if (!ehits)
0059   {
0060     std::cout << "PHG4CentralityReco::InitRun : cannot find G4HIT_EPD, will not use sEPD centrality" << std::endl;
0061   }
0062
0063   if (_centrality_calibration_params.exist_string_param("description"))
0064   {
0065     if (Verbosity() >= 1)
0066     {
0067       std::cout << "PHG4CentralityReco::InitRun : Centrality calibration description : " << std::endl
0068                 << "    ";
0069       std::cout << _centrality_calibration_params.get_string_param("description") << std::endl;
0070     }
0071     // search for possible centile definitions
0072     for (int n = 0; n < 101; n++)
0073     {
0074       std::ostringstream s1;
0075       s1 << "epd_centile_" << n;
0076       if (_centrality_calibration_params.exist_double_param(s1.str().c_str()))
0077       {
0078         _cent_cal_epd[_centrality_calibration_params.get_double_param(s1.str().c_str())] = n;
0079         if (Verbosity() >= 2)
0080         {
0081           std::cout << "PHG4CentralityReco::InitRun : sEPD centrality calibration, centile " << n << "% is " << _centrality_calibration_params.get_double_param(s1.str().c_str()) << std::endl;
0082         }
0083       }
0084     }
0085     for (int n = 0; n < 101; n++)
0086     {
0087       std::ostringstream s2;
0088       s2 << "mbd_centile_" << n;
0089       if (_centrality_calibration_params.exist_double_param(s2.str().c_str()))
0090       {
0091         _cent_cal_mbd[_centrality_calibration_params.get_double_param(s2.str().c_str())] = n;
0092         if (Verbosity() >= 2)
0093         {
0094           std::cout << "PHG4CentralityReco::InitRun : MBD centrality calibration, centile " << n << "% is " << _centrality_calibration_params.get_double_param(s2.str().c_str()) << std::endl;
0095         }
0096       }
0097     }
0098     for (int n = 0; n < 101; n++)
0099     {
0100       std::ostringstream s3;
0101       s3 << "bimp_centile_" << n;
0102       if (_centrality_calibration_params.exist_double_param(s3.str().c_str()))
0103       {
0104         _cent_cal_bimp[_centrality_calibration_params.get_double_param(s3.str().c_str())] = n;
0105         if (Verbosity() >= 2)
0106         {
0107           std::cout << "PHG4CentralityReco::InitRun : b (impact parameter) centrality calibration, centile " << n << "% is " << _centrality_calibration_params.get_double_param(s3.str().c_str()) << std::endl;
0108         }
0109       }
0110     }
0111   }
0112   else
0113   {
0114     std::cout << "PHG4CentralityReco::InitRun : no centrality calibration found!" << std::endl;
0115   }
0116
0117   return Fun4AllReturnCodes::EVENT_OK;
0118 }
0119
0120 int PHG4CentralityReco::process_event(PHCompositeNode *topNode)
0121 {
0122   _bimp = 101;
0123   auto event_header = findNode::getClass<EventHeaderv1>(topNode, "EventHeader");
0124   if (event_header)
0125   {
0126     _bimp = event_header->get_floatval("bimp");
0127
0128     if (Verbosity() >= 5)
0129     {
0130       std::cout << "PHG4CentralityReco::process_event : Hijing impact parameter b = " << _bimp << std::endl;
0131     }
0132   }
0133   else
0134   {
0135     if (Verbosity() >= 5)
0136     {
0137       std::cout << "PHG4CentralityReco::process_event : No Hijing impact parameter info, setting b = 101" << std::endl;
0138     }
0139   }
0140
0141   _mbd_N = 0;
0142   _mbd_S = 0;
0143   _mbd_NS = 0;
0144
0145   auto bhits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_BBC");
0146
0147   if (bhits)
0148   {
0149     auto brange = bhits->getHits();
0150     for (auto it = brange.first; it != brange.second; ++it)
0151     {
0152       if ((it->second->get_t(0) > -50) && (it->second->get_t(1) < 50))
0153       {
0154         _mbd_NS += it->second->get_edep();
0155         unsigned int id = it->second->get_layer();
0156         if ((id & 0x40U) == 0)
0157         {
0158           _mbd_N += it->second->get_edep();
0159         }
0160         else
0161         {
0162           _mbd_S += it->second->get_edep();
0163         }
0164       }
0165     }
0166
0167     if (Verbosity() >= 5)
0168     {
0169       std::cout << "PHG4CentralityReco::process_event : MBD Sum Charge N / S / N+S = " << _mbd_N << " / " << _mbd_S << " / " << _mbd_NS << std::endl;
0170     }
0171   }
0172   else
0173   {
0174     if (Verbosity() >= 5)
0175     {
0176       std::cout << "PHG4CentralityReco::process_event : No MBD info, setting all Sum Charges = 0" << std::endl;
0177     }
0178   }
0179
0180   _epd_N = 0;
0181   _epd_S = 0;
0182   _epd_NS = 0;
0183
0184   auto ehits = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_EPD");
0185
0186   if (ehits)
0187   {
0188     auto erange = ehits->getHits();
0189     for (auto it = erange.first; it != erange.second; ++it)
0190     {
0191       if ((it->second->get_t(0) > -50) && (it->second->get_t(1) < 50))
0192       {
0193         _epd_NS += it->second->get_edep();
0194
0195         int sim_tileid = (it->second->get_hit_id() >> PHG4HitDefs::hit_idbits);
0196         int this_arm = EPDDefs::get_arm(sim_tileid);
0197
0198         if (this_arm == 1)
0199         {
0200           _epd_N += it->second->get_edep();
0201         }
0202         else
0203         {
0204           _epd_S += it->second->get_edep();
0205         }
0206       }
0207     }
0208
0209     if (Verbosity() >= 5)
0210     {
0211       std::cout << "PHG4CentralityReco::process_event : sEPD Sum Energy N / S / N+S = " << _epd_N << " / " << _epd_S << " / " << _epd_NS << std::endl;
0212     }
0213   }
0214   else
0215   {
0216     if (Verbosity() >= 5)
0217     {
0218       std::cout << "PHG4CentralityReco::process_event : No sEPD info, setting all Sum Energies = 0" << std::endl;
0219     }
0220   }
0221
0222   if (Verbosity() >= 1)
0223   {
0224     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;
0225   }
0226
0227   if (_do_centrality_calibration)
0228   {
0229     // sEPD centrality
0230     float low_epd_val = -10000;
0231     float high_epd_val = 10000;
0232     int low_epd_centile = -1;
0233     int high_epd_centile = -1;
0234
0235     for (auto &it : _cent_cal_epd)
0236     {
0237       float signal = it.first;
0238       int cent = it.second;
0239
0240       if (signal < _epd_NS && signal > low_epd_val)
0241       {
0242         low_epd_val = signal;
0243         low_epd_centile = cent;
0244       }
0245       if (signal > _epd_NS && signal < high_epd_val)
0246       {
0247         high_epd_val = signal;
0248         high_epd_centile = cent;
0249       }
0250
0251     }  // close iterate through sEPD cuts
0252
0253     if (low_epd_centile >= 0 && high_epd_centile >= 0)
0254     {
0255       _epd_cent = (low_epd_centile + high_epd_centile) / 2.0;
0256       if (Verbosity() >= 10)
0257       {
0258         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;
0259       }
0260     }
0261     else
0262     {
0263       _epd_cent = 101;
0264       if (Verbosity() >= 5)
0265       {
0266         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;
0267       }
0268     }
0269
0270     // MBD centrality
0271     float low_mbd_val = -10000;
0272     float high_mbd_val = 10000;
0273     int low_mbd_centile = -1;
0274     int high_mbd_centile = -1;
0275
0276     for (auto &it : _cent_cal_mbd)
0277     {
0278       float signal = it.first;
0279       int cent = it.second;
0280
0281       if (signal < _mbd_NS && signal > low_mbd_val)
0282       {
0283         low_mbd_val = signal;
0284         low_mbd_centile = cent;
0285       }
0286       if (signal > _mbd_NS && signal < high_mbd_val)
0287       {
0288         high_mbd_val = signal;
0289         high_mbd_centile = cent;
0290       }
0291
0292     }  // close iterate through MBD cuts
0293
0294     if (low_mbd_centile >= 0 && high_mbd_centile >= 0)
0295     {
0296       _mbd_cent = (low_mbd_centile + high_mbd_centile) / 2.0;
0297       if (Verbosity() >= 10)
0298       {
0299         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;
0300       }
0301     }
0302     else
0303     {
0304       _mbd_cent = 101;
0305       if (Verbosity() >= 5)
0306       {
0307         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;
0308       }
0309     }
0310
0311     // b (impact parameter) centrality
0312     float low_bimp_val = -1;
0313     float high_bimp_val = 10000;
0314     int low_bimp_centile = -1;
0315     int high_bimp_centile = -1;
0316
0317     for (auto &it : _cent_cal_bimp)
0318     {
0319       float signal = it.first;
0320       int cent = it.second;
0321
0322       if (signal < _bimp && signal > low_bimp_val)
0323       {
0324         low_bimp_val = signal;
0325         low_bimp_centile = cent;
0326       }
0327       if (signal > _bimp && signal < high_bimp_val)
0328       {
0329         high_bimp_val = signal;
0330         high_bimp_centile = cent;
0331       }
0332
0333     }  // close iterate through bimp cuts
0334
0335     if (low_bimp_centile >= 0 && high_bimp_centile >= 0)
0336     {
0337       _bimp_cent = (low_bimp_centile + high_bimp_centile) / 2.0;
0338       if (Verbosity() >= 10)
0339       {
0340         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;
0341       }
0342     }
0343     else
0344     {
0345       _bimp_cent = 101;
0346       if (Verbosity() >= 5)
0347       {
0348         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;
0349       }
0350     }
0351
0352   }  // close centrality calibration
0353
0354   FillNode(topNode);
0355
0356   return Fun4AllReturnCodes::EVENT_OK;
0357 }
0358
0359 int PHG4CentralityReco::End(PHCompositeNode * /*topNode*/)
0360 {
0361   return Fun4AllReturnCodes::EVENT_OK;
0362 }
0363
0364 void PHG4CentralityReco::FillNode(PHCompositeNode *topNode)
0365 {
0366   CentralityInfo *cent = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0367   if (!cent)
0368   {
0369     std::cout << " ERROR -- can't find CentralityInfo node after it should have been created" << std::endl;
0370     return;
0371   }
0372   else
0373   {
0374     cent->set_quantity(CentralityInfo::PROP::mbd_N, _mbd_N);
0375     cent->set_quantity(CentralityInfo::PROP::mbd_S, _mbd_S);
0376     cent->set_quantity(CentralityInfo::PROP::mbd_NS, _mbd_NS);
0377     cent->set_quantity(CentralityInfo::PROP::epd_N, _epd_N);
0378     cent->set_quantity(CentralityInfo::PROP::epd_S, _epd_S);
0379     cent->set_quantity(CentralityInfo::PROP::epd_NS, _epd_NS);
0380     cent->set_quantity(CentralityInfo::PROP::bimp, _bimp);
0381
0382     cent->set_centile(CentralityInfo::PROP::epd_NS, _epd_cent);
0383     cent->set_centile(CentralityInfo::PROP::mbd_NS, _mbd_cent);
0384     cent->set_centile(CentralityInfo::PROP::bimp, _bimp_cent);
0385   }
0386 }
0387
0388 void PHG4CentralityReco::CreateNode(PHCompositeNode *topNode)
0389 {
0390   PHNodeIterator iter(topNode);
0391
0392   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0393   if (!dstNode)
0394   {
0395     std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0396     throw std::runtime_error("Failed to find DST node in PHG4CentralityReco::CreateNode");
0397   }
0398
0399   PHNodeIterator dstiter(dstNode);
0400   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "GLOBAL"));
0401   if (!DetNode)
0402   {
0403     DetNode = new PHCompositeNode("GLOBAL");
0404     dstNode->addNode(DetNode);
0405   }
0406
0407   CentralityInfo *cent = new CentralityInfov1();
0408
0409   PHIODataNode<PHObject> *centNode = new PHIODataNode<PHObject>(cent, "CentralityInfo", "PHObject");
0410   DetNode->addNode(centNode);
0411 }