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/PHG4TpcCylinderGeom.h>
0018 #include <g4detectors/PHG4TpcCylinderGeomContainer.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<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0050   if (!m_geom_container)
0051   {
0052     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << 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   TrkrHitSetContainer::ConstRange hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0165   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0166        hitsetitr != hitsetrange.second;
0167        ++hitsetitr)
0168   {
0169     TrkrHitSet *hitset = hitsetitr->second;
0170     int side = TpcDefs::getSide(hitsetitr->first);
0171
0172     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0173     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0174          hitr != hitrangei.second;
0175          ++hitr)
0176     {
0177       int it = TpcDefs::getTBin(hitr->first);
0178
0179       if (side == 0)
0180       {
0181         m_itHist_0->Fill(it);
0182       }
0183       else
0184       {
0185         m_itHist_1->Fill(it);
0186       }
0187     }
0188   }
0189
0190   double itMeanContent_0 = 0.0;
0191   double itMeanContent_1 = 0.0;
0192
0193   for (int i = 1; i <= m_time_samples_max; i++)
0194   {
0195     itMeanContent_0 += m_itHist_0->GetBinContent(i);
0196     itMeanContent_1 += m_itHist_1->GetBinContent(i);
0197   }
0198
0199   itMeanContent_0 = 1.0 * itMeanContent_0 / m_time_samples_max;
0200   itMeanContent_1 = 1.0 * itMeanContent_1 / m_time_samples_max;
0201
0202   m_itHist_0->GetXaxis()->SetRange(320, m_time_samples_max - 0.5);
0203   double itMax_0 = m_itHist_0->GetBinCenter(m_itHist_0->GetMaximumBin());
0204   double itMaxContent_0 = m_itHist_0->GetMaximum();
0205   m_itHist_0->GetXaxis()->SetRange(0, 0);
0206
0207   m_itHist_1->GetXaxis()->SetRange(320, m_time_samples_max - 0.5);
0208   double itMax_1 = m_itHist_1->GetBinCenter(m_itHist_1->GetMaximumBin());
0209   double itMaxContent_1 = m_itHist_1->GetMaximum();
0210   m_itHist_1->GetXaxis()->SetRange(0, 0);
0211
0212   auto f0 = std::make_unique<TF1>("f0", "gausn(0)");
0213   f0->SetParameters(itMaxContent_0, itMax_0, 1);
0214   f0->SetParLimits(1, itMax_0 - 2, itMax_0 + 2);
0215   m_itHist_0->Fit(f0.get(), "Bq0");
0216
0217   auto f1 = std::make_unique<TF1>("f1", "gausn(0)");
0218   f1->SetParameters(itMaxContent_1, itMax_1, 1);
0219   f1->SetParLimits(1, itMax_1 - 2, itMax_1 + 2);
0220   m_itHist_1->Fit(f1.get(), "Bq0");
0221
0222
0223   if ((itMaxContent_0 / itMeanContent_0 >= 7 && itMaxContent_0 > 1000) || (itMaxContent_1 / itMeanContent_1 >= 7 && itMaxContent_1 > 1000))
0224   {
0225     m_laserEventInfo->setIsLaserEvent(true);
0226     m_laserEventInfo->setPeakSample(false, (int) itMax_0);
0227     m_laserEventInfo->setPeakWidth(false, f0->GetParameter(2));
0228     m_laserEventInfo->setPeakSample(true, (int) itMax_1);
0229     m_laserEventInfo->setPeakWidth(true, f1->GetParameter(2));
0230
0231     isLaserEvent = true;
0232     peakSample0 = (int) itMax_0;
0233     peakSample1 = (int) itMax_1;
0234     peakWidth0 = f0->GetParameter(2);
0235     peakWidth1 = f1->GetParameter(2);
0236   }
0237   else
0238   {
0239     m_laserEventInfo->setIsLaserEvent(false);
0240
0241     isLaserEvent = false;
0242     peakSample0 = -999;
0243     peakSample1 = -999;
0244     peakWidth0 = -999;
0245     peakWidth1 = -999;
0246   }
0247
0248   if (m_debug)
0249   {
0250     m_hitTree->Fill();
0251   }
0252
0253
0254   return Fun4AllReturnCodes::EVENT_OK;
0255 }
0256
0257 int LaserEventIdentifier::ResetEvent(PHCompositeNode * /*topNode*/)
0258 {
0259   m_itHist_0->Reset();
0260   m_itHist_1->Reset();
0261
0262   return Fun4AllReturnCodes::EVENT_OK;
0263 }
0264
0265 int LaserEventIdentifier::End(PHCompositeNode * /*topNode*/)
0266 {
0267   if (m_debug)
0268   {
0269     m_debugFile->cd();
0270     m_hitTree->Write();
0271     m_debugFile->Close();
0272   }
0273
0274   return Fun4AllReturnCodes::EVENT_OK;
0275 }