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