mbdAnalysis/src/mbdAnalysis.cc

0001 //#include <bbc/BbcPmtContainerV1.h>
0002 //#include <bbc/BbcGeom.h>
0003 #include <mbd/MbdPmtContainer.h>
0004 #include <mbd/MbdGeom.h>
0005 #include <mbd/MbdPmtHit.h>
0006 #include <phool/getClass.h>
0007 #include <mbd/MbdPmtContainerV1.h>
0008 #include <mbd/MbdOut.h>
0009
0010 #include <phool/PHCompositeNode.h>
0011
0012 //Fun4All
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/Fun4AllServer.h>
0015 #include <fun4all/Fun4AllHistoManager.h>
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/getClass.h>
0018 #include <phool/phool.h>
0019 #include <ffaobjects/EventHeader.h>
0020
0021 //ROOT stuff
0022 #include <TH1F.h>
0023 #include <TH2F.h>
0024 #include <TH3F.h>
0025 #include <TFile.h>
0026 #include <TLorentzVector.h>
0027 #include <TTree.h>
0028
0029 #include "mbdAnalysis.h"
0030
0031 #include <fun4all/Fun4AllReturnCodes.h>
0032
0033 #include <phool/PHCompositeNode.h>
0034 // Centrality
0035 #include <centrality/CentralityInfo.h>
0036 // Minimum Bias
0037 #include <calotrigger/MinimumBiasInfo.h>
0038
0039
0040 #include <ffaobjects/EventHeader.h>
0041 #include <ffaobjects/EventHeaderv1.h>
0042
0043 //____________________________________________________________________________..
0044 mbdAnalysis::mbdAnalysis(const std::string &name):
0045  SubsysReco("poo")
0046     ,T(nullptr)
0047     ,Outfile(name)
0048 {
0049   std::cout << "mbdAnalysis::mbdAnalysis(const std::string &name) Calling ctor" << std::endl;
0050 }
0051
0052 //____________________________________________________________________________..
0053 mbdAnalysis::~mbdAnalysis()
0054 {
0055   std::cout << "mbdAnalysis::~mbdAnalysis() Calling dtor" << std::endl;
0056 }
0057
0058 //____________________________________________________________________________..
0059 int mbdAnalysis::Init(PHCompositeNode *topNode)
0060 {
0061
0062   out = new TFile(Outfile.c_str(),"RECREATE");
0063
0064   T = new TTree("T","T");
0065   T -> Branch("pmtcharge",&pmtcharge);
0066   T -> Branch("pmttime",&pmttime);
0067   T -> Branch("pmtx",&pmtx);
0068   T -> Branch("pmty",&pmty);
0069   T -> Branch("pmtz",&pmtz);
0070   T -> Branch("pmtr",&pmtr);
0071   T -> Branch("pmtphi",&pmtphi);
0072   T -> Branch("vertex",&vertex);
0073   T->Branch("cent",&cent);
0074   T -> Branch("evtPlaneAngle",&evtPlaneAngle);
0075
0076
0077   return Fun4AllReturnCodes::EVENT_OK;
0078 }
0079
0080 //____________________________________________________________________________..
0081 int mbdAnalysis::InitRun(PHCompositeNode *topNode)
0082 {
0083   std::cout << "mbdAnalysis::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0084   return Fun4AllReturnCodes::EVENT_OK;
0085 }
0086
0087 //____________________________________________________________________________..
0088
0089
0090
0091 int mbdAnalysis::process_event(PHCompositeNode *topNode)
0092 {
0093
0094
0095   //pmt information
0096   MbdPmtContainer *bbcpmts = findNode::getClass<MbdPmtContainer>(topNode,"MbdPmtContainer"); // bbc info
0097   if(!bbcpmts)
0098     {
0099       std::cout << "makeMBDTrees::process_event: Could not find MbdPmtContainer, aborting" << std::endl;
0100       return Fun4AllReturnCodes::ABORTEVENT;
0101     }
0102
0103   //MbdOut I think needed for vertex
0104   MbdOut *mbdout = findNode::getClass<MbdOut>(topNode,"MbdOut");
0105   if(!mbdout)
0106     {
0107       std::cout << "makeMBDTrees::process_event: Could not find MbdOut, aborting" << std::endl;
0108       return Fun4AllReturnCodes::ABORTEVENT;
0109     }
0110
0111   //Truth Informaiton
0112   EventHeader *evtHeader = findNode::getClass<EventHeader>(topNode,"EventHeader");
0113   if(!evtHeader)
0114     {
0115       std::cout << "makeMBDTrees::process_event: Could not find EventHeader, aborting" << std::endl;
0116       return Fun4AllReturnCodes::ABORTEVENT;
0117     }
0118 /*
0119   CentralityInfo *centInfo = findNode::getClass<CentralityInfo>(topNode,"CentralityInfo");
0120   if(!centInfo)
0121     {
0122       std::cout << PHWHERE << "caloTreeGen::process_event: Could not find node CentralityInfo" << std::endl;
0123       return Fun4AllReturnCodes::ABORTEVENT;
0124     }
0125   cent = centInfo->get_centile(CentralityInfo::PROP::mbd_NS);
0126 */
0127 /*
0128   MinimumBiasInfo *minBiasInfo = findNode::getClass<MinimumBiasInfo>(topNode,"MinimumBiasInfo");
0129   Bool_t isMinBias = (minBiasInfo) ? minBiasInfo->isAuAuMinimumBias() : false;
0130   if(!minBiasInfo || !isMinBias)
0131   {
0132     std::cout << PHWHERE << "caloTreeGen::process_event: is not MinimumBias" << std::endl;
0133     return Fun4AllReturnCodes::ABORTEVENT;
0134   }
0135   */
0136   //pmg geometry information
0137   MbdGeom *mbdgeom = findNode::getClass<MbdGeom>(topNode, "MbdGeom");
0138   if(!mbdgeom)
0139     {
0140       std::cout << "makeMBDTrees::process_event: Could not find MbdGeom, aborting" << std::endl;
0141       return Fun4AllReturnCodes::ABORTEVENT;
0142     }
0143
0144   //int nPMTs = bbcpmts -> get_npmt();        //size (should always be 128)
0145   int nPMTs = 128;        //size (should always be 128)
0146   vertex = mbdout->get_zvtx();
0147   std::cout << "NPMTS: "<< nPMTs << std::endl;
0148   evtPlaneAngle = evtHeader -> get_EventPlaneAngle();
0149   double m_cent = evtHeader->get_ImpactParameter();
0150   if(m_cent >= 0 && m_cent < 4.88) cent = 0.05;
0151   if(m_cent >= 4.88 && m_cent < 6.81) cent = 0.15;
0152   if(m_cent >= 6.81 && m_cent < 8.40) cent = 0.25;
0153   if(m_cent >= 8.40 && m_cent < 9.71) cent = 0.35;
0154   if(m_cent >= 9.71 && m_cent < 10.81) cent = 0.45;
0155   if(m_cent >= 10.81 && m_cent < 11.84) cent = 0.55;
0156   if(m_cent >= 11.84 && m_cent < 12.76) cent = 0.65;
0157   if(m_cent >= 12.76 && m_cent < 13.68) cent = 0.75;
0158   if(m_cent >= 13.68 && m_cent < 14.65) cent = 0.85;
0159
0160   for(int i = 0; i < nPMTs; i++)
0161     {
0162       MbdPmtHit* mbdpmt = bbcpmts -> get_pmt(i);        // grab ith pmt (changed from pmtID)
0163
0164       //store quantities of interest
0165       pmtcharge.push_back(mbdpmt -> get_q());
0166       pmttime.push_back(mbdpmt -> get_time());
0167       pmtx.push_back(mbdgeom->get_x(i));
0168       pmty.push_back(mbdgeom->get_y(i));
0169       pmtz.push_back(mbdgeom->get_z(i));
0170       pmtr.push_back(mbdgeom->get_r(i));
0171       pmtphi.push_back(mbdgeom->get_phi(i));
0172
0173     }
0174
0175   T -> Fill();
0176
0177
0178   return Fun4AllReturnCodes::EVENT_OK;
0179 }
0180
0181 //____________________________________________________________________________..
0182 int mbdAnalysis::ResetEvent(PHCompositeNode *topNode)
0183 {
0184   pmtcharge.clear();
0185   pmttime.clear();
0186   pmtx.clear();
0187   pmty.clear();
0188   pmtz.clear();
0189   pmtr.clear();
0190   pmtphi.clear();
0191   return Fun4AllReturnCodes::EVENT_OK;
0192 }
0193
0194 //____________________________________________________________________________..
0195 int mbdAnalysis::EndRun(const int runnumber)
0196 {
0197   std::cout << "mbdAnalysis::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0198   return Fun4AllReturnCodes::EVENT_OK;
0199 }
0200
0201 //____________________________________________________________________________..
0202 int mbdAnalysis::End(PHCompositeNode *topNode)
0203 {
0204   std::cout << "mbdAnalysis::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0205   out -> cd();
0206   T -> Write();
0207   out -> Close();
0208   return Fun4AllReturnCodes::EVENT_OK;
0209 }
0210
0211 //____________________________________________________________________________..
0212 int mbdAnalysis::Reset(PHCompositeNode *topNode)
0213 {
0214
0215   pmtcharge.clear();
0216   pmttime.clear();
0217   pmtx.clear();
0218   pmty.clear();
0219   pmtz.clear();
0220   pmtr.clear();
0221   pmtphi.clear();
0222
0223
0224   return Fun4AllReturnCodes::EVENT_OK;
0225 }
0226
0227 //____________________________________________________________________________..
0228 void mbdAnalysis::Print(const std::string &what) const
0229 {
0230   std::cout << "mbdAnalysis::Print(const std::string &what) const Printing info for " << what << std::endl;
0231 }