packages/mbd/MbdReco.cc

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