GlobalReco/Detinfo/Detinfo.cc

0001 #include "Detinfo.h"
0002
0003 #include <calobase/TowerInfoDefs.h>
0004 #include <epd/EPDDefs.h>
0005 #include <epd/EpdGeom.h>
0006
0007 #include <cdbobjects/CDBTTree.h> // for CDBTTree
0008 #include <ffamodules/CDBInterface.h>
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <fun4all/PHTFileServer.h>
0012
0013 #include <phool/getClass.h>
0014 #include <phool/PHCompositeNode.h>
0015
0016 #include <calobase/TowerInfoContainer.h>
0017 #include <calobase/TowerInfo.h>
0018
0019 #include <ffarawobjects/Gl1Packet.h>
0020 #include <phool/getClass.h>
0021
0022 #include <globalvertex/GlobalVertex.h>
0023 #include <globalvertex/GlobalVertexMap.h>
0024 #include <globalvertex/MbdVertex.h>
0025 #include <globalvertex/MbdVertexMapv1.h>
0026
0027 #include <calotrigger/MinimumBiasInfo.h>
0028 #include <centrality/CentralityInfo.h>
0029
0030 #include <eventplaneinfo/Eventplaneinfo.h>
0031 #include <eventplaneinfo/Eventplaneinfov1.h>
0032 #include <eventplaneinfo/EventplaneinfoMap.h>
0033
0034 #include <mbd/MbdGeom.h>
0035 #include <mbd/MbdOut.h>
0036 #include <mbd/MbdPmtContainer.h>
0037 #include <mbd/MbdPmtHit.h>
0038
0039 #include <zdcinfo/ZdcReco.h>
0040 #include <zdcinfo/Zdcinfo.h>
0041
0042 #include <boost/format.hpp>
0043
0044
0045 #include <TFile.h>
0046 #include <TH3F.h>
0047 #include <TH2F.h>
0048 #include <TTree.h>
0049 #include <TNtuple.h>
0050 #include <TProfile.h>
0051 #include <TVector3.h>
0052 #include <TH1.h>
0053
0054 #include <cassert>
0055 #include <fstream>
0056 #include <TGraph.h>
0057 #include <vector>
0058
0059
0060 #include <boost/format.hpp>
0061 #include <boost/math/special_functions/sign.hpp>
0062
0063
0064
0065 //____________________________________________________________________________..
0066 Detinfo::Detinfo(const std::string &name):
0067  SubsysReco(name)
0068 {
0069   _event = 0;
0070   _outfile_name = "qa.root";
0071
0072 }
0073
0074 //____________________________________________________________________________..
0075 Detinfo::~Detinfo()
0076 {
0077
0078 }
0079
0080 //____________________________________________________________________________..
0081 int Detinfo::Init(PHCompositeNode *topNode)
0082 {
0083
0084     std::cout << PHWHERE << " Opening file " << _outfile_name << std::endl;
0085
0086     PHTFileServer::get().open(_outfile_name, "RECREATE");
0087     PHTFileServer::get().cd(_outfile_name);
0088
0089     _event_tree = new TTree("event", "fwd => event info");
0090     _event_tree->Branch("event", &_event, "_event/I");
0091     _event_tree->Branch("zvertex", &thisvertex, "thisvertex/F");
0092     _event_tree->Branch("cent", &centrality_mbd_, "centrality_mbd_/F");
0093     _event_tree->Branch("psi2north", &psi2_north, "psi2_north/F");
0094     _event_tree->Branch("psi2south", &psi2_south, "psi2_south/F");
0095     _event_tree->Branch("mbdqsouth", &mbd_q_south, "mbd_q_south/F");
0096     _event_tree->Branch("mbdqnorth", &mbd_q_north, "mbd_q_north/F");
0097     _event_tree->Branch("sepdqsouth", &sepd_q_south, "sepd_q_south/F");
0098     _event_tree->Branch("sepdqnorth", &sepd_q_north, "sepd_q_north/F");
0099     _event_tree->Branch("zdcEsouth", &zdc_e_south, "zdc_e_south/F");
0100     _event_tree->Branch("zdcEnorth", &zdc_e_north, "zdc_e_north/F");
0101     _event_tree->Branch("is_min_bias", &is_min_bias);
0102     _event_tree->Branch("tile_e", &_t);
0103     _event_tree->Branch("tile_time", &_tt);
0104     _event_tree->Branch("pmt_e", &_p);
0105     _event_tree->Branch("pmt_time", &_pt);
0106     _event_tree->Branch("zdc_e", &_z);
0107     _event_tree->Branch("zdc_time", &_zt);
0108     h_triggerVec = new TH1D("h_triggerVec", "", 65, -0.5, 64.5);
0109
0110     return Fun4AllReturnCodes::EVENT_OK;
0111 }
0112
0113 //____________________________________________________________________________..
0114 int Detinfo::InitRun(PHCompositeNode *topNode)
0115 {
0116   //get sEPD tower mapping from CBD
0117   if (!m_overrideSEPDMapName) {
0118     m_sEPDMapName = "SEPD_CHANNELMAP";
0119   }
0120   if (!m_overrideSEPDFieldName) {
0121     m_sEPDfieldname = "epd_channel_map";
0122   }
0123   std::string calibdir = CDBInterface::instance()->getUrl(m_sEPDMapName);
0124   if (!calibdir.empty()) {
0125         cdbttree = new CDBTTree(calibdir);
0126  } else
0127  {
0128     std::cout << "Detinfo::::InitRun No SEPD mapping file for domain "<< m_sEPDMapName << " found" << std::endl;
0129     exit(1);
0130  }
0131
0132  v.clear();
0133  for (int i = 0; i < 768; i++)
0134  {
0135     int towerindex = cdbttree->GetIntValue(i, m_sEPDfieldname);
0136     if (towerindex == 999) //these are empty tiles
0137     {
0138         continue;
0139     }
0140     key = TowerInfoDefs::encode_epd(towerindex);
0141     v.push_back(key);
0142  }
0143   return Fun4AllReturnCodes::EVENT_OK;
0144 }
0145
0146 //____________________________________________________________________________..
0147 int Detinfo::process_event(PHCompositeNode *topNode)
0148 {
0149
0150     _event++;
0151
0152     _z.clear();
0153     _p.clear();
0154     _t.clear();
0155     _zt.clear();
0156     _pt.clear();
0157     _tt.clear();
0158
0159     mbd_q_south = 0.;
0160     mbd_q_north = 0.;
0161
0162     zdc_e_south = 0.;
0163     zdc_e_north = 0.;
0164
0165     sepd_q_south = 0.;
0166     sepd_q_north = 0.;
0167
0168     //get event plane map
0169     EventplaneinfoMap *epmap = findNode::getClass<EventplaneinfoMap>(topNode, "EventplaneinfoMap");
0170     if (!epmap)
0171     {
0172         std::cout << PHWHERE << "::ERROR - cannot find EventplaneinfoMap" << std::endl;
0173         exit(-1);
0174     }
0175
0176     if (epmap->empty())
0177     {
0178         return Fun4AllReturnCodes::ABORTEVENT;
0179     }
0180
0181     auto _EPDS = epmap->get(EventplaneinfoMap::sEPDS);
0182     auto _EPDN = epmap->get(EventplaneinfoMap::sEPDN);
0183
0184     //example psi_2 access
0185     std::pair<double, double> epdsouthQ2;
0186     std::pair<double, double> epdnorthQ2;
0187     epdsouthQ2 = _EPDS->get_qvector(2);
0188     epdnorthQ2 = _EPDN->get_qvector(2);
0189
0190     Qx_south = epdsouthQ2.first;
0191     Qy_south = epdsouthQ2.second;
0192     Qx_north = epdnorthQ2.first;
0193     Qy_north = epdnorthQ2.second;
0194
0195     Eventplaneinfo *epinfo = new Eventplaneinfov1();
0196     psi2_south = epinfo->GetPsi(Qx_south, Qy_south, 2);
0197     psi2_north = epinfo->GetPsi(Qx_north, Qy_north, 2);
0198
0199     CentralityInfo* m_CentInfo = nullptr;
0200     MinimumBiasInfo* _minimumbiasinfo = nullptr;
0201     m_CentInfo = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0202     if (!m_CentInfo)
0203     {
0204         std::cout << PHWHERE << "Error, can't find CentralityInfov1. No centrality info is filled" << std::endl;
0205          exit(1);
0206     }
0207
0208     if (m_CentInfo) centrality_mbd_ = m_CentInfo->get_centrality_bin(CentralityInfo::PROP::mbd_NS);
0209
0210     _minimumbiasinfo = findNode::getClass<MinimumBiasInfo>(topNode, "MinimumBiasInfo");
0211     if (!_minimumbiasinfo)
0212     {
0213         std::cout << "Error, can't find MinimumBiasInfo. No minimum bias info is filled" << std::endl;
0214         exit(1);
0215     }
0216
0217     is_min_bias = false;
0218     is_min_bias = (_minimumbiasinfo) ? _minimumbiasinfo->isAuAuMinimumBias() : false;
0219
0220
0221     MbdVertexMap *mbdmap =
0222           findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0223       MbdVertex *bvertex = nullptr;
0224       if (mbdmap)
0225       {
0226         for (MbdVertexMap::ConstIter mbditer = mbdmap->begin();
0227              mbditer != mbdmap->end(); ++mbditer)
0228         {
0229           bvertex = mbditer->second;
0230         }
0231         if (bvertex)
0232         {
0233             thisvertex = bvertex->get_z();
0234         }
0235       }
0236
0237     //--------------------------- trigger and GL1-------------------------------//
0238     Gl1Packet *gl1PacketInfo =
0239         findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0240     if (!gl1PacketInfo) {
0241       std::cout << PHWHERE << "Detinfo::process_event: GL1Packet node is missing"
0242                 << std::endl;
0243     }
0244
0245     uint64_t triggervec = 0;
0246     if (gl1PacketInfo) {
0247       triggervec = gl1PacketInfo->getScaledVector();
0248       for (int i = 0; i < 64; i++) {
0249         bool trig_decision = ((triggervec & 0x1U) == 0x1U);
0250
0251         if (trig_decision) {
0252           h_triggerVec->Fill(i);
0253         }
0254         triggervec = (triggervec >> 1U) & 0xffffffffU;
0255       }
0256       triggervec = gl1PacketInfo->getScaledVector();
0257     }
0258
0259
0260     TowerInfoContainer* towerinfosZDC;
0261     towerinfosZDC = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_ZDC");
0262     if(!towerinfosZDC)
0263     {
0264         std::cout << PHWHERE << ":: No TOWERINFO_CALIB_ZDC!" << std::endl; exit(1);
0265     }
0266
0267     MbdPmtHit *m_mbd_hit{nullptr};
0268     MbdPmtContainer *mbdpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0269      if (!mbdpmts)
0270      {
0271          std::cout << PHWHERE << "::ERROR - cannot find MbdPmtContainer" << std::endl;
0272          exit(1);
0273      }
0274
0275     TowerInfoContainer* towerinfosEPD = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_SEPD");
0276     if(!towerinfosEPD)
0277     {
0278         std::cout << PHWHERE << ":: No TOWERINFO_CALIB_SEPD!" << std::endl; exit(1);
0279     }
0280
0281     Zdcinfo *zdcinfo = findNode::getClass<Zdcinfo>(topNode, "Zdcinfo");
0282     if(!zdcinfo)
0283     {
0284         std::cout << PHWHERE << ":: No Zdcinfo!" << std::endl; exit(1);
0285     }
0286
0287     zdc_e_south = zdcinfo->get_zdc_energy(0);
0288     zdc_e_north = zdcinfo->get_zdc_energy(1);
0289
0290     MbdGeom *mbdgeom = findNode::getClass<MbdGeom>(topNode, "MbdGeom");
0291     if (!mbdgeom)
0292     {
0293        std::cout << PHWHERE << "::ERROR - cannot find MbdGeom" << std::endl;
0294         exit(-1);
0295     }
0296
0297
0298     if ((triggervec >> 0xAU) & 0x1U)
0299     {
0300         //mbd
0301         for (int i = 0; i < 128; i++)
0302         {
0303             m_mbd_hit = mbdpmts->get_pmt(i);
0304             _p.push_back(m_mbd_hit->get_q());
0305             _pt.push_back(m_mbd_hit->get_time());
0306
0307             //mbd charge sums
0308             int arm = mbdgeom->get_arm(i);
0309             if (arm == 0)
0310             {
0311                 mbd_q_south += m_mbd_hit->get_q();
0312             }
0313             else if (arm == 1)
0314             {
0315                 mbd_q_north += m_mbd_hit->get_q();
0316             }
0317         }
0318
0319         //sepd
0320         int nchannels_epd = towerinfosEPD->size();
0321
0322         float epd_e = 0.; float epd_t = 0.;
0323
0324         for (int channel = 0; channel < nchannels_epd;channel++)
0325         {
0326             epd_e = towerinfosEPD->get_tower_at_channel(channel)->get_energy();
0327             epd_t = towerinfosEPD->get_tower_at_channel(channel)->get_time_float();
0328             _t.push_back(epd_e);
0329             _tt.push_back(epd_t);
0330
0331             //sepd charge sums
0332             int arm = TowerInfoDefs::get_epd_arm(v[channel]);
0333             if (arm == 0)
0334             {
0335                 sepd_q_south += epd_e;
0336             }
0337             else if (arm == 1)
0338             {
0339                 sepd_q_north += epd_e;
0340             }
0341         }
0342
0343       //zdc
0344         int nchannels_zdc = towerinfosZDC->size();
0345         float zdc_t = 0.;
0346         float zdc_e = 0.;
0347
0348         for (int channel = 0; channel < nchannels_zdc;channel++)
0349         {
0350             zdc_t = towerinfosZDC->get_tower_at_channel(channel)->get_time_float();
0351             zdc_e = towerinfosZDC->get_tower_at_channel(channel)->get_energy();
0352             _z.push_back(zdc_e);
0353             _zt.push_back(zdc_t);
0354
0355         }
0356         _event_tree->Fill();
0357
0358     }//trigger
0359
0360
0361     std::cout<<_event<<std::endl;
0362
0363     return Fun4AllReturnCodes::EVENT_OK;
0364 }
0365
0366 //____________________________________________________________________________..
0367 int Detinfo::ResetEvent(PHCompositeNode *topNode)
0368 {
0369   return Fun4AllReturnCodes::EVENT_OK;
0370 }
0371
0372 //____________________________________________________________________________..
0373 int Detinfo::EndRun(const int runnumber)
0374 {
0375   return Fun4AllReturnCodes::EVENT_OK;
0376 }
0377
0378 //____________________________________________________________________________..
0379 int Detinfo::End(PHCompositeNode *topNode)
0380 {
0381     PHTFileServer::get().cd(_outfile_name);
0382     PHTFileServer::get().write(_outfile_name);
0383
0384     return Fun4AllReturnCodes::EVENT_OK;
0385
0386 }
0387
0388 //____________________________________________________________________________..
0389 int Detinfo::Reset(PHCompositeNode *topNode)
0390 {
0391   return Fun4AllReturnCodes::EVENT_OK;
0392 }