analysis/MBDinfo/MBDinfo.cc

0001
0002 // light weight code for MBD info
0003 // \author Ejiro Umaka
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 <mbd/MbdGeom.h>
0013 #include <mbd/MbdOut.h>
0014 #include <mbd/MbdPmtContainer.h>
0015 #include <mbd/MbdPmtHit.h>
0016
0017 #include <TFile.h>
0018 #include <TH3F.h>
0019 #include <TH2F.h>
0020 #include <TTree.h>
0021 #include <TNtuple.h>
0022
0023
0024 #include <cassert>
0025 #include <fstream>
0026 #include <TGraph.h>
0027 #include <vector>
0028 #include <TComplex.h>
0029
0030
0031 //____________________________________________________________________________..
0032 MBDinfo::MBDinfo(const std::string &name):
0033  SubsysReco(name)
0034 {
0035   _event = 0;
0036   _outfile_name = "MBDinfo.root";
0037
0038 }
0039
0040 //____________________________________________________________________________..
0041 MBDinfo::~MBDinfo()
0042 {
0043
0044 }
0045
0046 //____________________________________________________________________________..
0047 int MBDinfo::Init(PHCompositeNode *topNode)
0048 {
0049
0050     std::cout << PHWHERE << " Opening file " << _outfile_name << std::endl;
0051
0052     PHTFileServer::get().open(_outfile_name, "RECREATE");
0053     PHTFileServer::get().cd(_outfile_name);
0054     mbd = new TNtuple("mbd","mbd", "mbd_total_charge:mbd_south_total_charge:mbd_north_total_charge:mbd_south_raw_psi:mbd_north_raw_psi");
0055     return Fun4AllReturnCodes::EVENT_OK;
0056 }
0057
0058 //____________________________________________________________________________..
0059 int MBDinfo::InitRun(PHCompositeNode *topNode)
0060 {
0061   return Fun4AllReturnCodes::EVENT_OK;
0062 }
0063
0064 //____________________________________________________________________________..
0065 int MBDinfo::process_event(PHCompositeNode *topNode)
0066 {
0067
0068     _event++;
0069
0070     MbdPmtContainer *mbdpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0071      if (!mbdpmts)
0072      {
0073          std::cout << PHWHERE << "::ERROR - cannot find MbdPmtContainer" << std::endl;
0074          exit(-1);
0075      }
0076
0077      MbdGeom *mbdgeom = findNode::getClass<MbdGeom>(topNode, "MbdGeom");
0078      if (!mbdgeom)
0079      {
0080          std::cout << PHWHERE << "::ERROR - cannot find MbdGeom" << std::endl;
0081          exit(-1);
0082      }
0083
0084
0085     _f.clear();
0086
0087     float mbd_e_south = 0.;
0088     float mbd_e_north = 0.;
0089     float mbdQ = 0.;
0090     TComplex Q2_north(0.0,0.0);
0091     TComplex Q2_south(0.0,0.0);
0092
0093     //for (int ipmt = 0; ipmt < mbdpmts->get_npmt(); ipmt++) //issue with container size returning 0
0094     for (int ipmt = 0; ipmt < 128; ipmt++)
0095     {
0096         float mbd_q = mbdpmts->get_pmt(ipmt)->get_q();
0097         float phi = mbdgeom->get_phi(ipmt);
0098         int arm = mbdgeom->get_arm(ipmt);
0099         mbdQ += mbd_q;
0100
0101         if (arm == 0)
0102         {
0103             mbd_e_south += mbd_q;
0104             Q2_south += TComplex(mbd_q*cos(2*phi),mbd_q*sin(2*phi));
0105         }
0106
0107         else if (arm == 1)
0108         {
0109             mbd_e_north += mbd_q;
0110             Q2_north += TComplex(mbd_q*cos(2*phi),mbd_q*sin(2*phi));
0111
0112         }
0113     }
0114
0115     if((mbd_e_south > 0.0) && (mbd_e_north > 0.0))
0116     {
0117      _f.push_back(mbdQ);
0118      _f.push_back(mbd_e_south);
0119      _f.push_back(mbd_e_north);
0120      _f.push_back(atan2(Q2_south.Im(),Q2_south.Re())/2.0);
0121      _f.push_back(atan2(Q2_north.Im(),Q2_north.Re())/2.0);
0122       StoreMBDInfo(_f);
0123     }
0124
0125     std::cout<<_event<<std::endl;
0126     _f.clear();
0127
0128     return Fun4AllReturnCodes::EVENT_OK;
0129 }
0130
0131 //____________________________________________________________________________..
0132 int MBDinfo::ResetEvent(PHCompositeNode *topNode)
0133 {
0134   return Fun4AllReturnCodes::EVENT_OK;
0135 }
0136
0137 //____________________________________________________________________________..
0138 int MBDinfo::EndRun(const int runnumber)
0139 {
0140   return Fun4AllReturnCodes::EVENT_OK;
0141 }
0142
0143 //____________________________________________________________________________..
0144 int MBDinfo::End(PHCompositeNode *topNode)
0145 {
0146     PHTFileServer::get().cd(_outfile_name);
0147     PHTFileServer::get().write(_outfile_name);
0148
0149     return Fun4AllReturnCodes::EVENT_OK;
0150
0151 }
0152
0153 //____________________________________________________________________________..
0154 int MBDinfo::Reset(PHCompositeNode *topNode)
0155 {
0156   return Fun4AllReturnCodes::EVENT_OK;
0157 }
0158
0159
0160 double MBDinfo::StoreMBDInfo(std::vector < float > _m)
0161 {
0162
0163     for (int i = 0; i < 5; i++)
0164     {
0165         mtower_info[i] = _m[i];
0166     }
0167
0168
0169     mbd->Fill(mtower_info);
0170
0171     return 0;
0172
0173 }