packages/tpc/LaserEventIdentifier.cc

0001 #include "LaserEventIdentifier.h"
0002
0003 #include "LaserEventInfo.h"
0004 #include "LaserEventInfov2.h"
0005
0006 #include <trackbase/TpcDefs.h>
0007 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0008 #include <trackbase/TrkrHit.h>
0009 #include <trackbase/TrkrHitSet.h>
0010 #include <trackbase/TrkrHitSetContainer.h>
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/SubsysReco.h>  // for SubsysReco
0014
0015 #include <ffarawobjects/Gl1Packet.h>
0016
0017 #include <g4detectors/PHG4TpcGeom.h>
0018 #include <g4detectors/PHG4TpcGeomContainer.h>
0019
0020 #include <phool/PHCompositeNode.h>
0021 #include <phool/PHIODataNode.h>  // for PHIODataNode
0022 #include <phool/PHNode.h>        // for PHNode
0023 #include <phool/PHNodeIterator.h>
0024 #include <phool/PHObject.h>  // for PHObject
0025 #include <phool/getClass.h>
0026
0027 #include <Acts/Definitions/Units.hpp>
0028 #include <Acts/Surfaces/Surface.hpp>
0029
0030 #include <TF1.h>
0031 #include <TFile.h>
0032
0033 LaserEventIdentifier::LaserEventIdentifier(const std::string &name)
0034   : SubsysReco(name)
0035 {
0036 }
0037
0038 int LaserEventIdentifier::InitRun(PHCompositeNode *topNode)
0039 {
0040   // get node containing the digitized hits
0041   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0042   if (!m_hits)
0043   {
0044     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0045     return Fun4AllReturnCodes::ABORTRUN;
0046   }
0047
0048   m_geom_container =
0049       findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0050   if (!m_geom_container)
0051   {
0052     std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << std::endl;
0053     return Fun4AllReturnCodes::ABORTRUN;
0054   }
0055
0056   m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
0057                                                  "ActsGeometry");
0058   if (!m_tGeometry)
0059   {
0060     std::cout << PHWHERE
0061               << "ActsGeometry not found on node tree. Exiting"
0062               << std::endl;
0063     return Fun4AllReturnCodes::ABORTRUN;
0064   }
0065
0066   PHNodeIterator iter(topNode);
0067   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0068   if (!dstNode)
0069   {
0070     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0071     return Fun4AllReturnCodes::ABORTRUN;
0072   }
0073
0074   PHNodeIterator dstiter(dstNode);
0075   PHCompositeNode *DetNode =
0076       dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0077   if (!DetNode)
0078   {
0079     DetNode = new PHCompositeNode("TRKR");
0080     dstNode->addNode(DetNode);
0081   }
0082
0083   LaserEventInfo *laserEventInfo = new LaserEventInfov2();
0084
0085   PHIODataNode<PHObject> *laserEventInfoNode = new PHIODataNode<PHObject>(laserEventInfo, "LaserEventInfo", "PHObject");
0086   DetNode->addNode(laserEventInfoNode);
0087
0088   if (m_debug)
0089   {
0090     m_debugFile = new TFile(m_debugFileName.c_str(), "RECREATE");
0091   }
0092   float timeHistMax = m_time_samples_max;
0093   timeHistMax -= 0.5;
0094   m_itHist_0 = new TH1I("m_itHist_0", "side 0;it", m_time_samples_max, -0.5, timeHistMax);
0095   m_itHist_1 = new TH1I("m_itHist_1", "side 1;it", m_time_samples_max, -0.5, timeHistMax);
0096
0097   if (m_debug)
0098   {
0099     m_hitTree = new TTree("hitTree", "hitTree");
0100     m_hitTree->Branch("itHist_0", &m_itHist_0);
0101     m_hitTree->Branch("itHist_1", &m_itHist_1);
0102     m_hitTree->Branch("isLaserEvent", &isLaserEvent);
0103     m_hitTree->Branch("isGl1LaserEvent", &isGl1LaserEvent);
0104     m_hitTree->Branch("isGl1LaserPileupEvent", &isGl1LaserPileupEvent);
0105     m_hitTree->Branch("peakSample_0", &peakSample0);
0106     m_hitTree->Branch("peakSample_1", &peakSample1);
0107     m_hitTree->Branch("peakWidth_0", &peakWidth0);
0108     m_hitTree->Branch("peakWidth_1", &peakWidth1);
0109   }
0110
0111   return Fun4AllReturnCodes::EVENT_OK;
0112 }
0113
0114 int LaserEventIdentifier::process_event(PHCompositeNode *topNode)
0115 {
0116   m_laserEventInfo = findNode::getClass<LaserEventInfo>(topNode, "LaserEventInfo");
0117   if (!m_laserEventInfo)
0118   {
0119     std::cout << "no laser event info node" << std::endl;
0120     return Fun4AllReturnCodes::ABORTRUN;
0121   }
0122
0123   Gl1Packet *gl1pkt = findNode::getClass<Gl1Packet>(topNode, "GL1RAWHIT");
0124   if (!gl1pkt)
0125   {
0126     std::cout << "no GL1RAWHIT node" << std::endl;
0127     m_laserEventInfo->setIsGl1LaserEvent(false);
0128     m_laserEventInfo->setIsGl1LaserPileupEvent(false);
0129     //return Fun4AllReturnCodes::ABORTRUN;
0130   }
0131   else if(m_runnumber > 66153)
0132   {
0133     if ((gl1pkt->getGTMAllBusyVector() & (1<<14)) == 0)
0134     {
0135       m_laserEventInfo->setIsGl1LaserEvent(true);
0136       m_laserEventInfo->setIsGl1LaserPileupEvent(false);
0137       isGl1LaserEvent = true;
0138       isGl1LaserPileupEvent = false;
0139       prev_BCO = gl1pkt->getBCO();
0140     }
0141     else if ((gl1pkt->getBCO() - prev_BCO) < 350.0/30*16)
0142     {
0143       m_laserEventInfo->setIsGl1LaserEvent(false);
0144       m_laserEventInfo->setIsGl1LaserPileupEvent(true);
0145       isGl1LaserEvent = false;
0146       isGl1LaserPileupEvent = true;
0147       prev_BCO = 0;
0148     }
0149     else
0150     {
0151       m_laserEventInfo->setIsGl1LaserEvent(false);
0152       m_laserEventInfo->setIsGl1LaserPileupEvent(false);
0153       isGl1LaserEvent = false;
0154       isGl1LaserPileupEvent = false;
0155       prev_BCO = 0;
0156     }
0157   }
0158   else
0159   {
0160     m_laserEventInfo->setIsGl1LaserEvent(false);
0161     m_laserEventInfo->setIsGl1LaserPileupEvent(false);
0162   }
0163
0164
0165
0166   TrkrHitSetContainer::ConstRange hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0167   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0168        hitsetitr != hitsetrange.second;
0169        ++hitsetitr)
0170   {
0171     TrkrHitSet *hitset = hitsetitr->second;
0172     int side = TpcDefs::getSide(hitsetitr->first);
0173
0174     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0175     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0176          hitr != hitrangei.second;
0177          ++hitr)
0178     {
0179       int it = TpcDefs::getTBin(hitr->first);
0180
0181       if (side == 0)
0182       {
0183         m_itHist_0->Fill(it);
0184       }
0185       else
0186       {
0187         m_itHist_1->Fill(it);
0188       }
0189     }
0190   }
0191
0192   double itMeanContent_0 = 0.0;
0193   double itMeanContent_1 = 0.0;
0194
0195   for (int i = 1; i <= m_time_samples_max; i++)
0196   {
0197     itMeanContent_0 += m_itHist_0->GetBinContent(i);
0198     itMeanContent_1 += m_itHist_1->GetBinContent(i);
0199   }
0200
0201   itMeanContent_0 = 1.0 * itMeanContent_0 / m_time_samples_max;
0202   itMeanContent_1 = 1.0 * itMeanContent_1 / m_time_samples_max;
0203
0204   m_itHist_0->GetXaxis()->SetRange(320, m_time_samples_max - 0.5);
0205   double itMax_0 = m_itHist_0->GetBinCenter(m_itHist_0->GetMaximumBin());
0206   double itMaxContent_0 = m_itHist_0->GetMaximum();
0207   m_itHist_0->GetXaxis()->SetRange(0, 0);
0208
0209   m_itHist_1->GetXaxis()->SetRange(320, m_time_samples_max - 0.5);
0210   double itMax_1 = m_itHist_1->GetBinCenter(m_itHist_1->GetMaximumBin());
0211   double itMaxContent_1 = m_itHist_1->GetMaximum();
0212   m_itHist_1->GetXaxis()->SetRange(0, 0);
0213
0214   auto f0 = std::make_unique<TF1>("f0", "gausn(0)");
0215   f0->SetParameters(itMaxContent_0, itMax_0, 1);
0216   f0->SetParLimits(1, itMax_0 - 2, itMax_0 + 2);
0217   m_itHist_0->Fit(f0.get(), "Bq0");
0218
0219   auto f1 = std::make_unique<TF1>("f1", "gausn(0)");
0220   f1->SetParameters(itMaxContent_1, itMax_1, 1);
0221   f1->SetParLimits(1, itMax_1 - 2, itMax_1 + 2);
0222   m_itHist_1->Fit(f1.get(), "Bq0");
0223
0224
0225   if ((itMaxContent_0 / itMeanContent_0 >= 7 && itMaxContent_0 > 1000) || (itMaxContent_1 / itMeanContent_1 >= 7 && itMaxContent_1 > 1000))
0226   {
0227     m_laserEventInfo->setIsLaserEvent(true);
0228     m_laserEventInfo->setPeakSample(false, (int) itMax_0);
0229     m_laserEventInfo->setPeakWidth(false, f0->GetParameter(2));
0230     m_laserEventInfo->setPeakSample(true, (int) itMax_1);
0231     m_laserEventInfo->setPeakWidth(true, f1->GetParameter(2));
0232
0233     isLaserEvent = true;
0234     peakSample0 = (int) itMax_0;
0235     peakSample1 = (int) itMax_1;
0236     peakWidth0 = f0->GetParameter(2);
0237     peakWidth1 = f1->GetParameter(2);
0238   }
0239   else
0240   {
0241     m_laserEventInfo->setIsLaserEvent(false);
0242
0243     isLaserEvent = false;
0244     peakSample0 = -999;
0245     peakSample1 = -999;
0246     peakWidth0 = -999;
0247     peakWidth1 = -999;
0248   }
0249
0250   if (m_debug)
0251   {
0252     m_hitTree->Fill();
0253   }
0254
0255
0256   return Fun4AllReturnCodes::EVENT_OK;
0257 }
0258
0259 int LaserEventIdentifier::ResetEvent(PHCompositeNode * /*topNode*/)
0260 {
0261   m_itHist_0->Reset();
0262   m_itHist_1->Reset();
0263
0264   return Fun4AllReturnCodes::EVENT_OK;
0265 }
0266
0267 int LaserEventIdentifier::End(PHCompositeNode * /*topNode*/)
0268 {
0269   if (m_debug)
0270   {
0271     m_debugFile->cd();
0272     m_hitTree->Write();
0273     m_debugFile->Close();
0274   }
0275
0276   return Fun4AllReturnCodes::EVENT_OK;
0277 }