packages/mbd/MbdReco.cc

0001 #include "MbdReco.h"
0002 #include "MbdEvent.h"
0003 #include "MbdGeomV1.h"
0004 #include "MbdOutV2.h"
0005 #include "MbdPmtContainerV1.h"
0006 #include "MbdPmtSimContainerV1.h"
0007
0008 #include <globalvertex/MbdVertexMapv1.h>
0009 #include <globalvertex/MbdVertexv2.h>
0010
0011 #include <ffarawobjects/CaloPacket.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015
0016 #include <Event/Event.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHIODataNode.h>
0020 #include <phool/PHNode.h>
0021 #include <phool/PHNodeIterator.h>
0022 #include <phool/PHObject.h>
0023 #include <phool/PHRandomSeed.h>
0024 #include <phool/getClass.h>
0025 #include <phool/phool.h>
0026
0027 #include <ffaobjects/EventHeader.h>
0028 #include <ffarawobjects/CaloPacketContainer.h>
0029 #include <ffarawobjects/Gl1Packet.h>
0030
0031 #include <TF1.h>
0032
0033 //____________________________________________________________________________..
0034 MbdReco::MbdReco(const std::string &name)
0035   : SubsysReco(name)
0036 {
0037 }
0038
0039 //____________________________________________________________________________..
0040 int MbdReco::Init(PHCompositeNode * /*topNode*/)
0041 {
0042   m_gaussian = std::make_unique<TF1>("gaussian", "gaus", 0, 20);
0043   m_gaussian->FixParameter(2, m_tres);
0044
0045   m_mbdevent = std::make_unique<MbdEvent>(_calpass,_always_process_charge);
0046   if ( Verbosity()>0 )
0047   {
0048     m_mbdevent->Verbosity( Verbosity() );
0049   }
0050
0051   return Fun4AllReturnCodes::EVENT_OK;
0052 }
0053
0054 //____________________________________________________________________________..
0055 int MbdReco::InitRun(PHCompositeNode *topNode)
0056 {
0057   if (createNodes(topNode) == Fun4AllReturnCodes::ABORTEVENT)
0058   {
0059     return Fun4AllReturnCodes::ABORTEVENT;
0060   }
0061
0062   int ret = getNodes(topNode);
0063
0064   m_mbdevent->SetSim(_simflag);
0065   m_mbdevent->InitRun();
0066
0067   return ret;
0068 }
0069
0070 //____________________________________________________________________________..
0071 int MbdReco::process_event(PHCompositeNode *topNode)
0072 {
0073   getNodes(topNode);
0074
0075   if ( (m_mbdevent==nullptr && m_mbdraw==nullptr && m_mbdpacket[0]==nullptr && m_mbdpacket[1]==nullptr) || m_mbdpmts==nullptr )
0076   {
0077     static int counter = 0;
0078     if ( counter<2 )
0079     {
0080       std::cout << PHWHERE << " ERROR, didn't find mbdevent, mbdraw, or mbdpmts" << std::endl;
0081       counter++;
0082     }
0083     return Fun4AllReturnCodes::ABORTEVENT;  // missing an essential object in BBC/MBD
0084   }
0085
0086   if (m_mbdpacket[0] && m_mbdpacket[1] &&
0087       (m_mbdpacket[0]->getIdentifier() != 1001 || m_mbdpacket[1]->getIdentifier() != 1002))
0088   {
0089     static int counter = 0;
0090     if (counter < 100)
0091     {
0092       std::cout << PHWHERE << "packet 1001 and/or packet 1002 missing, bailing out" << std::endl;
0093       counter++;
0094     }
0095     return Fun4AllReturnCodes::EVENT_OK; // no mbd packets here
0096   }
0097
0098   // Process raw waveforms from real data
0099   if ( m_mbdevent!=nullptr || m_mbdraw!=nullptr || m_mbdpacket[0]==nullptr || m_mbdpacket[1]==nullptr)
0100   {
0101     int status = Fun4AllReturnCodes::EVENT_OK;
0102     if ( m_evtheader!=nullptr )
0103     {
0104       m_mbdevent->set_EventNumber( m_evtheader->get_EvtSequence() );
0105     }
0106     if ( m_event!=nullptr )
0107     {
0108       status = m_mbdevent->SetRawData(m_event, m_mbdpmts);
0109     }
0110     else if ( m_mbdraw!=nullptr || m_mbdpacket[0]!=nullptr || m_mbdpacket[1]!=nullptr)
0111     {
0112       if (m_mbdraw)
0113       {
0114     m_mbdpacket[0] = m_mbdraw->getPacketbyId(1001);
0115     m_mbdpacket[1] = m_mbdraw->getPacketbyId(1002);
0116       }
0117       status = m_mbdevent->SetRawData(m_mbdpacket, m_mbdpmts,m_gl1raw);
0118     }
0119
0120     if (status == Fun4AllReturnCodes::DISCARDEVENT )
0121     {
0122       static int counter = 0;
0123       if ( counter<3 )
0124       {
0125         std::cout << PHWHERE << " Warning, MBD discarding event " << std::endl;
0126         counter++;
0127       }
0128       return Fun4AllReturnCodes::DISCARDEVENT;
0129     }
0130     if (status == Fun4AllReturnCodes::ABORTEVENT )
0131     {
0132       static int counter = 0;
0133       if ( counter<3 )
0134       {
0135         std::cout << PHWHERE << " Warning, MBD aborting event " << std::endl;
0136         counter++;
0137       }
0138       return Fun4AllReturnCodes::ABORTEVENT;
0139     }
0140     if ( status == -1001 )
0141     {
0142       // calculating sampmax on this event
0143       return Fun4AllReturnCodes::DISCARDEVENT;
0144     }
0145     if (status < 0)
0146     {
0147       return Fun4AllReturnCodes::EVENT_OK;
0148     }
0149   }
0150
0151   // Calibrate from UNCALDST or recalibrate from DST
0152   if ( _calpass==3 )
0153   {
0154     m_mbdevent->ProcessRawPackets( m_mbdpmts );
0155   }
0156
0157   m_mbdevent->Calculate(m_mbdpmts, m_mbdout, topNode);
0158
0159   // For multiple global vertex
0160   if (m_mbdevent->get_bbcn(0) > 0 && m_mbdevent->get_bbcn(1) > 0 && _calpass==0 )
0161   {
0162     auto *vertex = new MbdVertexv2();
0163     vertex->set_t(m_mbdevent->get_bbct0());
0164     vertex->set_z(m_mbdevent->get_bbcz());
0165     vertex->set_z_err(0.6);
0166     vertex->set_t_err(m_tres);
0167     vertex->set_beam_crossing(0);
0168
0169     m_mbdvtxmap->insert(vertex);
0170
0171     // copy to globalvertex
0172   }
0173
0174   if (Verbosity() > 0)
0175   {
0176     std::cout << "mbd vertex z and t0 " << m_mbdevent->get_bbcz() << ", " << m_mbdevent->get_bbct0() << std::endl;
0177   }
0178
0179   return Fun4AllReturnCodes::EVENT_OK;
0180 }
0181
0182 //____________________________________________________________________________..
0183 int MbdReco::End(PHCompositeNode * /*unused*/)
0184 {
0185   m_mbdevent->End();
0186
0187   return Fun4AllReturnCodes::EVENT_OK;
0188 }
0189
0190 int MbdReco::createNodes(PHCompositeNode *topNode)
0191 {
0192   PHNodeIterator iter(topNode);
0193   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0194   if (!dstNode)
0195   {
0196     std::cout << PHWHERE << "DST Node missing doing nothing" << std::endl;
0197     return Fun4AllReturnCodes::ABORTRUN;
0198   }
0199
0200   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0201   if (!runNode)
0202   {
0203     std::cout << PHWHERE << "RUN Node missing doing nothing" << std::endl;
0204     return Fun4AllReturnCodes::ABORTRUN;
0205   }
0206
0207   PHNodeIterator dstiter(dstNode);
0208   PHCompositeNode *bbcNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "MBD"));
0209   if (!bbcNode)
0210   {
0211     bbcNode = new PHCompositeNode("MBD");
0212     dstNode->addNode(bbcNode);
0213   }
0214
0215   PHNodeIterator runiter(runNode);
0216   PHCompositeNode *bbcRunNode = dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode", "MBD"));
0217   if (!bbcRunNode)
0218   {
0219     bbcRunNode = new PHCompositeNode("MBD");
0220     runNode->addNode(bbcRunNode);
0221   }
0222
0223   m_mbdout = findNode::getClass<MbdOut>(bbcNode, "MbdOut");
0224   if (!m_mbdout)
0225   {
0226     m_mbdout = new MbdOutV2();
0227     PHIODataNode<PHObject> *MbdOutNode = new PHIODataNode<PHObject>(m_mbdout, "MbdOut", "PHObject");
0228     bbcNode->addNode(MbdOutNode);
0229   }
0230
0231   m_mbdpmts = findNode::getClass<MbdPmtSimContainerV1>(bbcNode, "MbdPmtContainer");
0232   if (!m_mbdpmts)
0233   {
0234     m_mbdpmts = new MbdPmtContainerV1();
0235
0236     PHIODataNode<PHObject> *MbdPmtContainerNode = new PHIODataNode<PHObject>(m_mbdpmts, "MbdPmtContainer", "PHObject");
0237     bbcNode->addNode(MbdPmtContainerNode);
0238   }
0239
0240   PHCompositeNode *globalNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "GLOBAL"));
0241   if (!globalNode)
0242   {
0243     globalNode = new PHCompositeNode("GLOBAL");
0244     dstNode->addNode(globalNode);
0245   }
0246
0247   m_mbdvtxmap = findNode::getClass<MbdVertexMap>(globalNode, "MbdVertexMap");
0248   if (!m_mbdvtxmap)
0249   {
0250     m_mbdvtxmap = new MbdVertexMapv1();
0251     PHIODataNode<PHObject> *VertexMapNode = new PHIODataNode<PHObject>(m_mbdvtxmap, "MbdVertexMap", "PHObject");
0252     globalNode->addNode(VertexMapNode);
0253   }
0254
0255   m_mbdgeom = findNode::getClass<MbdGeom>(runNode, "MbdGeom");
0256   if (!m_mbdgeom)
0257   {
0258     m_mbdgeom = new MbdGeomV1();
0259     PHIODataNode<PHObject> *MbdGeomNode = new PHIODataNode<PHObject>(m_mbdgeom, "MbdGeom", "PHObject");
0260     bbcRunNode->addNode(MbdGeomNode);
0261   }
0262
0263   return Fun4AllReturnCodes::EVENT_OK;
0264 }
0265
0266 int MbdReco::getNodes(PHCompositeNode *topNode)
0267 {
0268   // Get the bbc prdf data to mpcRawContent
0269   m_event = findNode::getClass<Event>(topNode, "PRDF");
0270   // std::cout << "event addr " << (unsigned int)m_event << endl;
0271
0272   // Get the raw data from event combined DST
0273   m_mbdraw = findNode::getClass<CaloPacketContainer>(topNode, "MBDPackets");
0274
0275   m_mbdpacket[0] = findNode::getClass<CaloPacket>(topNode,1001);
0276   m_mbdpacket[1] = findNode::getClass<CaloPacket>(topNode,1002);
0277   if (!m_event && !m_mbdraw && !m_mbdpacket[0] && !m_mbdpacket[1])
0278   {
0279     // not PRDF and not event combined DST, so we assume this is a sim file
0280     _simflag = 1;
0281
0282     static int counter = 0;
0283     if (counter < 1)
0284     {
0285       std::cout << PHWHERE << "Unable to get PRDF or Event Combined DST, assuming this is simulation" << std::endl;
0286       counter++;
0287     }
0288   }
0289
0290   // Get the raw gl1 data from event combined DST
0291   m_gl1raw = findNode::getClass<Gl1Packet>(topNode,14001);
0292   if (!m_gl1raw)
0293   {
0294     m_gl1raw = findNode::getClass<Gl1Packet>(topNode, "GL1Packet");
0295   }
0296   /*
0297   if ( !m_gl1raw )
0298   {
0299     cout << PHWHERE << " Gl1Packet node not found on node tree" << endl;
0300   }
0301   */
0302
0303
0304   // MbdPmtContainer
0305   m_mbdpmts = findNode::getClass<MbdPmtContainer>(topNode, "MbdPmtContainer");
0306   if (!m_mbdpmts)
0307   {
0308     std::cout << PHWHERE << " MbdPmtContainer node not found on node tree" << std::endl;
0309     return Fun4AllReturnCodes::ABORTEVENT;
0310   }
0311
0312   m_mbdvtxmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0313   if (!m_mbdvtxmap)
0314   {
0315     std::cout << PHWHERE << "MbdVertexMap node not found on node tree" << std::endl;
0316     return Fun4AllReturnCodes::ABORTEVENT;
0317   }
0318
0319   m_evtheader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0320   if (!m_evtheader )
0321   {
0322     static int ctr = 0;
0323     if ( ctr<4 )
0324     {
0325       std::cout << PHWHERE << " EventHeader node not found on node tree" << std::endl;
0326       ctr++;
0327     }
0328   }
0329
0330   return Fun4AllReturnCodes::EVENT_OK;
0331 }