analysis/EventplaneValidation/MBDinfo.cc

0001 //Author Ejiro Umaka
0002 //sample code for extracting MBD eventplane info
0003
0004 #include "MBDinfo.h"
0005
0006 #include <fun4all/Fun4AllReturnCodes.h>
0007 #include <fun4all/PHTFileServer.h>
0008
0009 #include <phool/getClass.h>
0010 #include <phool/PHCompositeNode.h>
0011
0012 #include <eventplaneinfo/Eventplaneinfo.h>
0013 #include <eventplaneinfo/EventplaneinfoMap.h>
0014
0015 #include <globalvertex/GlobalVertex.h>
0016 #include <globalvertex/GlobalVertexMap.h>
0017
0018 #include <centrality/CentralityInfo.h>
0019 #include <centrality/CentralityInfov1.h>
0020
0021
0022 #include <mbd/MbdGeom.h>
0023 #include <mbd/MbdOut.h>
0024 #include <mbd/MbdPmtContainer.h>
0025 #include <mbd/MbdPmtHit.h>
0026
0027 #include <TFile.h>
0028 #include <TH3F.h>
0029 #include <TH2F.h>
0030 #include <TTree.h>
0031 #include <TNtuple.h>
0032
0033
0034 #include <cassert>
0035 #include <fstream>
0036 #include <TGraph.h>
0037 #include <vector>
0038
0039
0040
0041
0042 //____________________________________________________________________________..
0043 MBDinfo::MBDinfo(const std::string &name):
0044  SubsysReco(name)
0045 {
0046   _event = 0;
0047   _outfile_name = "MBDinfo.root";
0048
0049 }
0050
0051 //____________________________________________________________________________..
0052 MBDinfo::~MBDinfo()
0053 {
0054
0055 }
0056
0057 //____________________________________________________________________________..
0058 int MBDinfo::Init(PHCompositeNode *topNode)
0059 {
0060
0061     std::cout << PHWHERE << " Opening file " << _outfile_name << std::endl;
0062
0063     PHTFileServer::get().open(_outfile_name, "RECREATE");
0064     PHTFileServer::get().cd(_outfile_name);
0065     mbd = new TNtuple("mbd","mbd", "vertex:cent:tQ:sQ:nQ:spsi1:spsi1_cor:spsi2:spsi2_cor:spsi3:spsi3_cor:npsi1:npsi1_cor:npsi2:npsi2_cor:npsi3:npsi3_cor");
0066     return Fun4AllReturnCodes::EVENT_OK;
0067 }
0068
0069 //____________________________________________________________________________..
0070 int MBDinfo::InitRun(PHCompositeNode *topNode)
0071 {
0072   return Fun4AllReturnCodes::EVENT_OK;
0073 }
0074
0075 //____________________________________________________________________________..
0076 int MBDinfo::process_event(PHCompositeNode *topNode)
0077 {
0078
0079     _event++;
0080     _f.clear();
0081
0082      MbdPmtContainer *mbdpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0083      if (!mbdpmts)
0084      {
0085          std::cout << PHWHERE << "::ERROR - cannot find MbdPmtContainer" << std::endl;
0086          exit(-1);
0087      }
0088
0089      MbdGeom *mbdgeom = findNode::getClass<MbdGeom>(topNode, "MbdGeom");
0090      if (!mbdgeom)
0091      {
0092          std::cout << PHWHERE << "::ERROR - cannot find MbdGeom" << std::endl;
0093          exit(-1);
0094      }
0095
0096
0097     GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0098     if (!vertexmap)
0099     {
0100       std::cout << PHWHERE << " Fatal Error - GlobalVertexMap node is missing" << std::endl;
0101       return Fun4AllReturnCodes::ABORTRUN;
0102     }
0103
0104    if (!(vertexmap->empty()))
0105     {
0106         GlobalVertex *vtx = vertexmap->begin()->second;
0107         if (vtx)
0108         {
0109            thisvertex = vtx->get_z();
0110         }
0111     }
0112
0113   _f.push_back(thisvertex);
0114
0115
0116 if(_dosim){
0117       CentralityInfov1 *cent = findNode::getClass<CentralityInfov1>(topNode, "CentralityInfo");
0118       if (!cent)
0119       {
0120           std::cout << "can't find CentralityInfo node " << std::endl;
0121           exit(1);
0122
0123       }
0124
0125     float cent_index = cent->get_centile(CentralityInfo::PROP::bimp);
0126     _f.push_back(cent_index);
0127 }
0128 else
0129 {
0130   _f.push_back(-990.0);
0131 }
0132
0133
0134     float mbd_e_south = 0.;
0135     float mbd_e_north = 0.;
0136     float mbdQ = 0.;
0137
0138
0139     for (int ipmt = 0; ipmt < 128; ipmt++)
0140     {
0141         float mbd_q = mbdpmts->get_pmt(ipmt)->get_q();
0142         int arm = mbdgeom->get_arm(ipmt);
0143         mbdQ += mbd_q;
0144
0145         if (arm == 0)
0146         {
0147             mbd_e_south += mbd_q;
0148         }
0149
0150         else if (arm == 1)
0151         {
0152             mbd_e_north += mbd_q;
0153
0154         }
0155     }
0156
0157
0158
0159     double sep1 = NAN;
0160     double nep1 = NAN;
0161     double sep_shift1 = NAN;
0162     double nep_shift1 = NAN;
0163
0164
0165     double sep2 = NAN;
0166     double nep2 = NAN;
0167     double sep_shift2 = NAN;
0168     double nep_shift2 = NAN;
0169
0170
0171     double sep3 = NAN;
0172     double nep3 = NAN;
0173     double sep_shift3 = NAN;
0174     double nep_shift3 = NAN;
0175
0176
0177     EventplaneinfoMap *epmap = findNode::getClass<EventplaneinfoMap>(topNode, "EventplaneinfoMap");
0178     if (!epmap)
0179     {
0180        std::cout << PHWHERE << "::ERROR - cannot find EventplaneinfoMap" << std::endl;
0181          exit(-1);
0182     }
0183
0184     if (epmap->empty())
0185     {
0186 //      std::cout << "Empty epmap - this event falls outside abs(60cm) or has no determined vertex" << std::endl;
0187       return Fun4AllReturnCodes::ABORTEVENT;
0188     }
0189
0190     auto mbds = epmap->get(EventplaneinfoMap::MBDS);
0191
0192     sep1 = mbds->get_psi(1);
0193     sep2 = mbds->get_psi(2);
0194     sep3 = mbds->get_psi(3);
0195     sep_shift1 = mbds->get_shifted_psi(1);
0196     sep_shift2 = mbds->get_shifted_psi(2);
0197     sep_shift3 = mbds->get_shifted_psi(3);
0198
0199     auto mbdn = epmap->get(EventplaneinfoMap::MBDN);
0200
0201     nep1 = mbdn->get_psi(1);
0202     nep2 = mbdn->get_psi(2);
0203     nep3 = mbdn->get_psi(3);
0204     nep_shift1 = mbdn->get_shifted_psi(1);
0205     nep_shift2 = mbdn->get_shifted_psi(2);
0206     nep_shift3 = mbdn->get_shifted_psi(3);
0207
0208      _f.push_back(mbdQ);
0209      _f.push_back(mbd_e_south);
0210      _f.push_back(mbd_e_north);
0211
0212     _f.push_back(sep1);
0213     _f.push_back(sep_shift1);
0214
0215     _f.push_back(sep2);
0216     _f.push_back(sep_shift2);
0217
0218     _f.push_back(sep3);
0219     _f.push_back(sep_shift3);
0220
0221
0222     _f.push_back(nep1);
0223     _f.push_back(nep_shift1);
0224
0225     _f.push_back(nep2);
0226     _f.push_back(nep_shift2);
0227
0228     _f.push_back(nep3);
0229     _f.push_back(nep_shift3);
0230
0231     StoreMBDInfo(_f);
0232
0233     std::cout<<_event<<std::endl;
0234
0235      //QVector info
0236
0237     //std::pair<double, double> Qsouth;
0238     //Qsouth = mbds->get_qvector(2);// second order qvector in south
0239
0240 //    std::pair<double, double> Qnorth;
0241 //    Qnorth = mbdn->get_qvector(2);// second order qvector in north
0242
0243     //float qxsouth = Qsouth.first;
0244     //float qysouth = Qsouth.second;
0245
0246     _f.clear();
0247
0248     return Fun4AllReturnCodes::EVENT_OK;
0249 }
0250
0251 //____________________________________________________________________________..
0252 int MBDinfo::ResetEvent(PHCompositeNode *topNode)
0253 {
0254   return Fun4AllReturnCodes::EVENT_OK;
0255 }
0256
0257 //____________________________________________________________________________..
0258 int MBDinfo::EndRun(const int runnumber)
0259 {
0260   return Fun4AllReturnCodes::EVENT_OK;
0261 }
0262
0263 //____________________________________________________________________________..
0264 int MBDinfo::End(PHCompositeNode *topNode)
0265 {
0266     PHTFileServer::get().cd(_outfile_name);
0267     PHTFileServer::get().write(_outfile_name);
0268
0269     return Fun4AllReturnCodes::EVENT_OK;
0270
0271 }
0272
0273 //____________________________________________________________________________..
0274 int MBDinfo::Reset(PHCompositeNode *topNode)
0275 {
0276   return Fun4AllReturnCodes::EVENT_OK;
0277 }
0278
0279
0280 double MBDinfo::StoreMBDInfo(std::vector < float > _m)
0281 {
0282
0283     for (int i = 0; i < 17; i++)
0284     {
0285         mtower_info[i] = _m[i];
0286     }
0287
0288
0289     mbd->Fill(mtower_info);
0290
0291     return 0;
0292
0293 }