centralmembraneanalysis/TpcCentralMembranePlotter/TpcCentralMembranePlotter.cc

0001 #include "TpcCentralMembranePlotter.h"
0002
0003 #include <tpc/LaserEventInfo.h>
0004
0005
0006 #include <trackbase/CMFlashDifferenceContainerv1.h>
0007 #include <trackbase/CMFlashDifferencev1.h>
0008 #include <trackbase/LaserClusterContainerv1.h>
0009 #include <trackbase/LaserClusterv3.h>
0010 #include <trackbase/TpcDefs.h>
0011 #include <trackbase/TrkrHitTruthAssoc.h>
0012 #include <trackbase/TrkrHitSetContainer.h>
0013 #include <trackbase/TrkrHitSet.h>
0014 #include <trackbase/TrkrHit.h>
0015
0016 #include <g4main/PHG4Hit.h>
0017 #include <g4main/PHG4HitContainer.h>
0018
0019 #include <g4detectors/PHG4TpcGeom.h>
0020 #include <g4detectors/PHG4TpcGeomContainer.h>
0021
0022 #include <ffaobjects/EventHeader.h>
0023
0024 #include <cdbobjects/CDBTTree.h>
0025
0026 #include <fun4all/DBInterface.h>
0027 #include <fun4all/Fun4AllReturnCodes.h>
0028
0029 #include <phool/PHCompositeNode.h>
0030 #include <phool/getClass.h>
0031 #include <phool/phool.h>
0032
0033 #include <odbc++/resultset.h>
0034 #include <odbc++/statement.h>
0035
0036 #include <TCanvas.h>
0037 #include <TF1.h>
0038 #include <TFile.h>
0039 #include <TGraph.h>
0040 #include <TH2.h>
0041 #include <TH3.h>
0042 #include <TLine.h>
0043 #include <TLegend.h>
0044 #include <TPaveText.h>
0045 #include <TStyle.h>
0046 #include <TTree.h>
0047 #include <TVector3.h>
0048
0049 #include <boost/format.hpp>
0050 #include <boost/math/special_functions/lambert_w.hpp>
0051
0052
0053 #include <cmath>
0054 #include <format>
0055 #include <iomanip>
0056 #include <set>
0057 #include <string>
0058 #include <vector>
0059
0060 //________________________________
0061 TpcCentralMembranePlotter::TpcCentralMembranePlotter(const std::string &name)
0062   : SubsysReco(name)
0063 {
0064 }
0065
0066 //____________________________________________
0067 int TpcCentralMembranePlotter::InitRun(PHCompositeNode *topNode)
0068 {
0069   for (int s = 0; s < 2; s++)
0070   {
0071     float zmin = (s == 1 ? -40 : 30);
0072     float zmax = (s == 1 ? -30 : 40);
0073
0074     m_hits_hardware_layervsphibin[s] = new TH2D((boost::format("m_hits_hardware_layervsphibin_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;#phi bin;Layer") %(s == 1 ? "North" : "South")).str().c_str(), 2600, 0, 2600, 50, 6, 56);
0075     m_hits_hardware_timebinvsphibin[s] = new TH2D((boost::format("m_hits_hardware_timebinvsphibin_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;#phi bin;Time bin") %(s == 1 ? "North" : "South")).str().c_str(), 2600, 0, 2600, 40,310,350);
0076     m_hits_hardware_timebinvslayer[s] = new TH2D((boost::format("m_hits_hardware_timebinvslayer_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;Layer;Time bin") %(s == 1 ? "North" : "South")).str().c_str(), 50, 6, 56, 40,310,350);
0077
0078     m_clusters_r_vs_phi[s] = new TH2D((boost::format("m_clusters_r_vs_phi_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;#phi;R [cm]") %(s == 1 ? "North" : "South")).str().c_str(), 6280, 0, 6.28, 200, 30, 80);
0079     m_clusters_z_vs_phi[s] = new TH2D((boost::format("m_clusters_z_vs_phi_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;#phi;Z [cm]") %(s == 1 ? "North" : "South")).str().c_str(), 6280, 0, 6.28, 200, zmin, zmax);
0080     m_clusters_z_vs_r[s] = new TH2D((boost::format("m_clusters_z_vs_r_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;R [cm];Z [cm]") %(s == 1 ? "North" : "South")).str().c_str(), 1000, 30, 80, 1000, zmin, zmax);
0081
0082     m_hits_r_vs_phi[s] = new TH2D((boost::format("m_hits_r_vs_phi_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;#phi;R [cm]") %(s == 1 ? "North" : "South")).str().c_str(), 6280, 0, 6.28, 200, 30, 80);
0083     m_hits_z_vs_phi[s] = new TH2D((boost::format("m_hits_z_vs_phi_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;#phi;Z [cm]") %(s == 1 ? "North" : "South")).str().c_str(), 6280, 0, 6.28, 1000, zmin, zmax);
0084     m_hits_z_vs_r[s] = new TH2D((boost::format("m_hits_z_vs_r_%s") %(s == 1 ? "North" : "South")).str().c_str(), (boost::format("TPC %s;R [cm];Z [cm]") %(s == 1 ? "North" : "South")).str().c_str(), 1000, 30, 80, 1000, zmin, zmax);
0085
0086
0087   }
0088
0089
0090   int ret = GetNodes(topNode);
0091   return ret;
0092 }
0093
0094 //__________________________________________
0095 int TpcCentralMembranePlotter::process_event(PHCompositeNode *topNode)
0096 {
0097   eventHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0098   if (!eventHeader)
0099   {
0100     std::cout << PHWHERE << " EventHeader Node missing, abort" << std::endl;
0101     return Fun4AllReturnCodes::ABORTRUN;
0102   }
0103
0104   m_runnumber = eventHeader->get_RunNumber();
0105
0106   if (m_useHeader && eventHeader->get_EvtSequence() == 0)
0107   {
0108     m_useHeader = false;
0109   }
0110
0111   if (m_useHeader)
0112   {
0113     m_event_index = eventHeader->get_EvtSequence();
0114   }
0115
0116   if (!m_correctedCMcluster_map || m_correctedCMcluster_map->size() < 1000)
0117   {
0118     if (!m_useHeader)
0119     {
0120       m_event_index++;
0121     }
0122     return Fun4AllReturnCodes::EVENT_OK;
0123   }
0124
0125   if ((eventHeader->get_RunNumber() > 66153 && !m_laserEventInfo->isGl1LaserEvent()) || (eventHeader->get_RunNumber() <= 66153 && !m_laserEventInfo->isLaserEvent()))
0126   {
0127     std::cout << "Not a laser event, skipping " << m_event_index << std::endl;
0128     return Fun4AllReturnCodes::EVENT_OK;
0129   }
0130
0131   m_nClusters += m_correctedCMcluster_map->size();
0132   m_nEvents++;
0133
0134   std::cout << "TpcCentralMembranePlotter::process_event - Processing event " << m_event_index << std::endl;
0135   std::cout << m_nClusters << std::endl;
0136
0137   auto clusrange = m_correctedCMcluster_map->getClusters();
0138
0139   for (auto cmitr = clusrange.first; cmitr != clusrange.second; ++cmitr)
0140   {
0141     const auto &[cmkey, cmclus_orig] = *cmitr;
0142     LaserCluster *cmclus = cmclus_orig;
0143     const unsigned int clusadc = cmclus->getAdc();
0144     bool side = (bool) TpcDefs::getSide(cmkey);
0145
0146     if (cmclus->getNLayers() < m_nLayerCut)
0147     {
0148       continue;
0149     }
0150
0151     // ******* Fill histogram of hits ******* //
0152     int nHitsInClus = cmclus->getNhits();
0153     for(int i=0; i< nHitsInClus; ++i)
0154     {
0155         const LaserClusterHitInfo hit= cmclus->getHit(i);
0156
0157         TrkrHitTruthAssoc::MMap g4hit_map;
0158         m_hitTruthAssoc->getG4Hits(hit.hitsetkey, hit.hitkey, g4hit_map);
0159
0160         // const bool has_truth = !g4hit_map.empty();
0161
0162         // std::cout
0163         // << (has_truth ? "SIGNAL " : "NOISE  ")
0164         // << "hitsetkey " << hit.hitsetkey
0165         // << " trkrid " << static_cast<int>(TrkrDefs::getTrkrId(hit.hitsetkey))
0166         // << " layer " << static_cast<int>(TrkrDefs::getLayer(hit.hitsetkey))
0167         // << " side " << static_cast<int>(TpcDefs::getSide(hit.hitsetkey))
0168         // << " pad " << static_cast<int>(TpcDefs::getPad(hit.hitkey))
0169         // << " tbin " << static_cast<int>(TpcDefs::getTBin(hit.hitkey))
0170         // << " adc " << hit.adc
0171         // << " ng4 " << g4hit_map.size()
0172         // << std::endl;
0173
0174
0175         int hitside = TpcDefs::getSide(hit.hitsetkey);
0176
0177         double hitadc = hit.adc;
0178         int hitlayer = TrkrDefs::getLayer(hit.hitsetkey);
0179         int hitiphi = TpcDefs::getPad(hit.hitkey);
0180         int hitit = TpcDefs::getTBin(hit.hitkey);
0181
0182         PHG4TpcGeom *layer_geom = m_geom_container->GetLayerCellGeom(hitlayer);
0183         double hitenvradius = layer_geom->get_radius();
0184         double hitenvphi = layer_geom->get_phi(hitiphi, hitside);
0185         double hitenv_x = hitenvradius * cos(hitenvphi);
0186         double hitenv_y = hitenvradius * sin(hitenvphi);
0187         double hitenv_z = 0.0;
0188         Acts::Vector3 hitenv(hitenv_x, hitenv_y, hitenv_z);
0189         Acts::Vector3 global = m_tGeometry->transformTpcEnvelopeToWorld(hitenv);
0190         double hitglob_x = global[0];
0191         double hitglob_y = global[1];
0192         double hitglob_z = global[2];
0193         double hitglob_r = sqrt(pow(hitglob_x,2) + pow(hitglob_y,2));
0194         double hitglob_phi = atan2(hitglob_y,hitglob_x);
0195
0196         if (hitenvphi < 0.0)
0197         {
0198           hitenvphi += 2 * M_PI;
0199         }
0200         if (hitglob_phi < 0.0)
0201         {
0202           hitglob_phi += 2 * M_PI;
0203         }
0204
0205         m_hits_hardware_timebinvsphibin[hitside]->Fill(hitiphi, hitit, hitadc);
0206         m_hits_hardware_timebinvslayer[hitside]->Fill(hitlayer, hitit, hitadc);
0207         m_hits_hardware_layervsphibin[hitside]->Fill(hitiphi, hitlayer, hitadc);
0208
0209         m_hits_z_vs_phi[hitside]->Fill(hitglob_phi, hitglob_z, hitadc);
0210         m_hits_z_vs_r[hitside]->Fill(hitglob_r, hitglob_z, hitadc);
0211         m_hits_r_vs_phi[hitside]->Fill(hitglob_phi, hitglob_r, hitadc);
0212
0213
0214     }
0215
0216     Acts::Vector3 pos = m_laserClusterHelper.getClusterCentroid(cmclus);
0217     if (m_dcc_in_module_edge)
0218     {
0219       pos = m_distortionCorrection.get_corrected_position(pos, m_dcc_in_module_edge);
0220     }
0221     if (m_dcc_in_static)
0222     {
0223       pos = m_distortionCorrection.get_corrected_position(pos, m_dcc_in_static);
0224     }
0225
0226     TVector3 tmp_pos(pos[0], pos[1], pos[2]);
0227
0228     // if(cmclus->getSDWeightedLayer() > 0.5)
0229     // {
0230     //   continue;
0231     // }
0232
0233     double phi2pimod = tmp_pos.Phi();
0234     if (phi2pimod < 0.0)
0235     {
0236       phi2pimod += 2 * M_PI;
0237     }
0238     m_clusters_r_vs_phi[side]->Fill(phi2pimod, tmp_pos.Perp(), clusadc);
0239     m_clusters_z_vs_phi[side]->Fill(phi2pimod, tmp_pos.Z(), clusadc);
0240     m_clusters_z_vs_r[side]->Fill(tmp_pos.Perp(), tmp_pos.Z(), clusadc);
0241
0242   }
0243
0244   return Fun4AllReturnCodes::EVENT_OK;
0245 }
0246
0247 int TpcCentralMembranePlotter::End(PHCompositeNode * /*topNode*/)
0248 {
0249   TFile *outputfile = new TFile(m_outputfile.c_str(), "RECREATE");
0250   outputfile->cd();
0251   for (int s = 0; s < 2; s++)
0252   {
0253     m_hits_hardware_layervsphibin[s]->Write();
0254     m_hits_hardware_timebinvsphibin[s]->Write();
0255     m_hits_hardware_timebinvslayer[s]->Write();
0256
0257   }
0258
0259   for (int s = 0; s < 2; s++)
0260   {
0261     m_hits_r_vs_phi[s]->Write();
0262     m_hits_z_vs_phi[s]->Write();
0263     m_hits_z_vs_r[s]->Write();
0264   }
0265
0266   for (int s = 0; s < 2; s++)
0267   {
0268     m_clusters_r_vs_phi[s]->Write();
0269     m_clusters_z_vs_phi[s]->Write();
0270     m_clusters_z_vs_r[s]->Write();
0271   }
0272
0273   outputfile->Close();
0274
0275   return Fun4AllReturnCodes::EVENT_OK;
0276 }
0277
0278 //______________________________________
0279 int TpcCentralMembranePlotter::GetNodes(PHCompositeNode *topNode)
0280 {
0281   m_correctedCMcluster_map = findNode::getClass<LaserClusterContainer>(topNode, "LASER_CLUSTER");
0282   if (!m_correctedCMcluster_map)
0283   {
0284     std::cout << PHWHERE << "CORRECTED_CM_CLUSTER Node missing, abort." << std::endl;
0285     return Fun4AllReturnCodes::ABORTRUN;
0286   }
0287
0288   m_dcc_in_module_edge = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerModuleEdge");
0289   if (m_dcc_in_module_edge)
0290   {
0291     std::cout << "TpcCentralMembranePlotter::GetNodes - found TPC distortion correction container module edge" << std::endl;
0292   }
0293
0294   m_dcc_in_static = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, "TpcDistortionCorrectionContainerStatic");
0295   if (m_dcc_in_static)
0296   {
0297     std::cout << "TpcCentralMembranePlotter::GetNodes - found TPC distortion correction container static" << std::endl;
0298   }
0299
0300   PHNodeIterator iter(topNode);
0301
0302   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0303   if (!dstNode)
0304   {
0305     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0306     return Fun4AllReturnCodes::ABORTRUN;
0307   }
0308
0309   auto *flashDiffContainer = findNode::getClass<CMFlashDifferenceContainerv1>(topNode, "CM_FLASH_DIFFERENCES");
0310   if (!flashDiffContainer)
0311   {
0312     PHNodeIterator dstIter(dstNode);
0313     PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "TRKR"));
0314     if (!DetNode)
0315     {
0316       DetNode = new PHCompositeNode("TRKR");
0317       dstNode->addNode(DetNode);
0318     }
0319
0320     flashDiffContainer = new CMFlashDifferenceContainerv1;
0321     PHIODataNode<PHObject> *CMFlashDifferenceNode = new PHIODataNode<PHObject>(flashDiffContainer, "CM_FLASH_DIFFERENCES", "PHObject");
0322     DetNode->addNode(CMFlashDifferenceNode);
0323   }
0324
0325   const std::string dcc_out_node_name = "TpcDistortionCorrectionContainerAverage";
0326   m_dcc_out = findNode::getClass<TpcDistortionCorrectionContainer>(topNode, dcc_out_node_name);
0327   if (!m_dcc_out)
0328   {
0329     auto *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0330     if (!runNode)
0331     {
0332       std::cout << "TpcCentralMembranePlotter::InitRun = RUN Node missing, quitting" << std::endl;
0333       return Fun4AllReturnCodes::ABORTRUN;
0334     }
0335
0336     std::cout << "TpcCentralMembranePlotter::GetNodes - creating TpcDistortionCorrectionContainer in node " << dcc_out_node_name << std::endl;
0337     m_dcc_out = new TpcDistortionCorrectionContainer;
0338     auto *node = new PHDataNode<TpcDistortionCorrectionContainer>(m_dcc_out, dcc_out_node_name);
0339     runNode->addNode(node);
0340   }
0341
0342   m_tGeometry = findNode::getClass<ActsGeometry>(topNode,"ActsGeometry");
0343     if(!m_tGeometry)
0344     {
0345       std::cout << PHWHERE << "ERROR: Can't find node ActsGeometry" << std::endl;
0346     }
0347
0348   m_geom_container = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0349   if (!m_geom_container)
0350   {
0351     std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << std::endl;
0352     return Fun4AllReturnCodes::ABORTRUN;
0353   }
0354
0355   m_laserEventInfo = findNode::getClass<LaserEventInfo>(topNode, "LaserEventInfo");
0356   if (!m_laserEventInfo)
0357   {
0358     std::cout << PHWHERE << "ERROR: Can't find node LaserEventInfo" << std::endl;
0359     return Fun4AllReturnCodes::ABORTRUN;
0360   }
0361
0362   m_hitTruthAssoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0363   if (!m_hitTruthAssoc)
0364   {
0365     std::cout << PHWHERE << "ERROR: Can't find TRKR_HITTRUTHASSOC" << std::endl;
0366     return Fun4AllReturnCodes::ABORTRUN;
0367   }
0368
0369   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0370   if (!m_hits)
0371   {
0372     std::cout << PHWHERE << "ERROR: Can't find TRKR_HITSET" << std::endl;
0373     return Fun4AllReturnCodes::ABORTRUN;
0374   }
0375
0376   m_laserClusterHelper.set_useZ(0);
0377   m_laserClusterHelper.loadNodes(topNode);
0378
0379   return Fun4AllReturnCodes::EVENT_OK;
0380 }
0381
0382 //_____________________________
0383 void TpcCentralMembranePlotter::setOutputfile(const std::string &outputfile)
0384 {
0385   m_outputfile = outputfile;
0386 }