QA/Tracking/TpcSiliconQA.cc

0001 #include "TpcSiliconQA.h"
0002
0003 #include <trackbase/TpcDefs.h>
0004 #include <trackbase/TrkrDefs.h>
0005
0006 #include <trackbase_historic/TrackSeed.h>
0007 #include <trackbase_historic/TrackSeedContainer.h>
0008 #include <trackbase_historic/TrackSeedHelper.h>
0009
0010 #include <qautils/QAHistManagerDef.h>
0011 #include <qautils/QAUtil.h>
0012
0013 #include <fun4all/Fun4AllHistoManager.h>
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015 #include <fun4all/SubsysReco.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/getClass.h>
0019
0020 #include <TH1.h>
0021 #include <TH2.h>
0022
0023 #include <boost/format.hpp>
0024 #include <iomanip>
0025
0026 //____________________________________________________________________________..
0027 TpcSiliconQA::TpcSiliconQA(const std::string& name)
0028   : SubsysReco(name)
0029 {
0030 }
0031
0032 //____________________________________________________________________________..
0033 int TpcSiliconQA::InitRun(PHCompositeNode* /*topNode*/)
0034 {
0035   createHistos();
0036
0037   return Fun4AllReturnCodes::EVENT_OK;
0038 }
0039
0040 //____________________________________________________________________________..
0041 int TpcSiliconQA::process_event(PHCompositeNode* topNode)
0042 {
0043   auto* hm = QAHistManagerDef::getHistoManager();
0044   assert(hm);
0045
0046   auto* silseedmap = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0047   if (!silseedmap)
0048   {
0049     std::cout << "Silicon seed map not found, aborting event" << std::endl;
0050     return Fun4AllReturnCodes::ABORTEVENT;
0051   }
0052   auto* tpcseedmap = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0053   if (!tpcseedmap)
0054   {
0055     std::cout << "TPC seed map not found, aborting event" << std::endl;
0056     return Fun4AllReturnCodes::ABORTEVENT;
0057   }
0058
0059   for (const auto& silseed : *silseedmap)
0060   {
0061     if (!silseed)
0062     {
0063       continue;
0064     }
0065
0066     m_crossing = (float) silseed->get_crossing();
0067     h_crossing->Fill(m_crossing);
0068
0069     {
0070       const auto position = TrackSeedHelper::get_xyz(silseed);
0071       m_silseedx = position.x();
0072       m_silseedy = position.y();
0073       m_silseedz = position.z();
0074       m_silseedphi = silseed->get_phi();
0075       m_silseedeta = silseed->get_eta();
0076     }
0077
0078     for (const auto& tpcseed : *tpcseedmap)
0079     {
0080       if (!tpcseed)
0081       {
0082         continue;
0083       }
0084
0085       const auto position = TrackSeedHelper::get_xyz(tpcseed);
0086       m_tpcseedx = position.x();
0087       m_tpcseedy = position.y();
0088       m_tpcseedz = position.z();
0089       m_tpcseedphi = tpcseed->get_phi();
0090       m_tpcseedeta = tpcseed->get_eta();
0091
0092       h_phiDiff[0]->Fill(m_tpcseedphi - m_silseedphi);
0093       h_etaDiff[0]->Fill(m_tpcseedeta - m_silseedeta);
0094       h_xDiff[0]->Fill(m_tpcseedx - m_silseedx);
0095       h_yDiff[0]->Fill(m_tpcseedy - m_silseedy);
0096       h_zDiff[0]->Fill(m_tpcseedz - m_silseedz);
0097
0098       if (m_tpcseedeta > 0 && m_silseedeta > 0)
0099       {
0100         h_phiDiff[4]->Fill(m_tpcseedphi - m_silseedphi);
0101         h_etaDiff[4]->Fill(m_tpcseedeta - m_silseedeta);
0102         h_xDiff[4]->Fill(m_tpcseedx - m_silseedx);
0103         h_yDiff[4]->Fill(m_tpcseedy - m_silseedy);
0104         h_zDiff[4]->Fill(m_tpcseedz - m_silseedz);
0105       }
0106       else if (m_tpcseedeta < 0 && m_silseedeta < 0)
0107       {
0108         h_phiDiff[5]->Fill(m_tpcseedphi - m_silseedphi);
0109         h_etaDiff[5]->Fill(m_tpcseedeta - m_silseedeta);
0110         h_xDiff[5]->Fill(m_tpcseedx - m_silseedx);
0111         h_yDiff[5]->Fill(m_tpcseedy - m_silseedy);
0112         h_zDiff[5]->Fill(m_tpcseedz - m_silseedz);
0113       }
0114
0115       if (abs(m_tpcseedx - m_silseedx) > m_xcut || abs(m_tpcseedy - m_silseedy) > m_ycut)
0116       {
0117         continue;
0118       }
0119
0120       h_phiDiff[1]->Fill(m_tpcseedphi - m_silseedphi);
0121       h_etaDiff[1]->Fill(m_tpcseedeta - m_silseedeta);
0122       h_xDiff[1]->Fill(m_tpcseedx - m_silseedx);
0123       h_yDiff[1]->Fill(m_tpcseedy - m_silseedy);
0124       h_zDiff[1]->Fill(m_tpcseedz - m_silseedz);
0125
0126       if (abs(m_tpcseedeta - m_silseedeta) > m_etacut)
0127       {
0128         continue;
0129       }
0130
0131       h_phiDiff[2]->Fill(m_tpcseedphi - m_silseedphi);
0132       h_etaDiff[2]->Fill(m_tpcseedeta - m_silseedeta);
0133       h_xDiff[2]->Fill(m_tpcseedx - m_silseedx);
0134       h_yDiff[2]->Fill(m_tpcseedy - m_silseedy);
0135       h_zDiff[2]->Fill(m_tpcseedz - m_silseedz);
0136
0137       if (abs(m_tpcseedphi - m_silseedphi) > m_phicut)
0138       {
0139         continue;
0140       }
0141
0142       h_phiDiff[3]->Fill(m_tpcseedphi - m_silseedphi);
0143       h_etaDiff[3]->Fill(m_tpcseedeta - m_silseedeta);
0144       h_xDiff[3]->Fill(m_tpcseedx - m_silseedx);
0145       h_yDiff[3]->Fill(m_tpcseedy - m_silseedy);
0146       h_zDiff[3]->Fill(m_tpcseedz - m_silseedz);
0147
0148       if (m_tpcseedeta > 0 && m_silseedeta > 0)
0149       {
0150         h_phiDiff[6]->Fill(m_tpcseedphi - m_silseedphi);
0151         h_etaDiff[6]->Fill(m_tpcseedeta - m_silseedeta);
0152         h_xDiff[6]->Fill(m_tpcseedx - m_silseedx);
0153         h_yDiff[6]->Fill(m_tpcseedy - m_silseedy);
0154         h_zDiff[6]->Fill(m_tpcseedz - m_silseedz);
0155       }
0156       else if (m_tpcseedeta < 0 && m_silseedeta < 0)
0157       {
0158         h_phiDiff[7]->Fill(m_tpcseedphi - m_silseedphi);
0159         h_etaDiff[7]->Fill(m_tpcseedeta - m_silseedeta);
0160         h_xDiff[7]->Fill(m_tpcseedx - m_silseedx);
0161         h_yDiff[7]->Fill(m_tpcseedy - m_silseedy);
0162         h_zDiff[7]->Fill(m_tpcseedz - m_silseedz);
0163       }
0164     }
0165   }
0166
0167   m_event++;
0168   return Fun4AllReturnCodes::EVENT_OK;
0169 }
0170
0171 int TpcSiliconQA::EndRun(const int /*runnumber*/)
0172 {
0173   return Fun4AllReturnCodes::EVENT_OK;
0174 }
0175
0176 std::string TpcSiliconQA::getHistoPrefix() const
0177 {
0178   return std::string("h_") + Name() + std::string("_");
0179 }
0180
0181 void TpcSiliconQA::createHistos()
0182 {
0183   auto* hm = QAHistManagerDef::getHistoManager();
0184   assert(hm);
0185
0186   std::stringstream stream1;
0187   std::stringstream stream2;
0188   std::stringstream stream3;
0189   std::stringstream stream4;
0190   stream1 << std::fixed << std::setprecision(2) << m_xcut;
0191   stream2 << std::fixed << std::setprecision(2) << m_ycut;
0192   stream3 << std::fixed << std::setprecision(2) << m_etacut;
0193   stream4 << std::fixed << std::setprecision(2) << m_phicut;
0194
0195   std::vector<std::string> cutNames = {"", "_xyCut", "_etaCut", "_phiCut", "North", "South", "NorthAllCuts", "SouthAllCuts"};
0196   std::vector<std::string> cutVals = {"All Track Seeds",
0197                                       std::string("|xdiff| < " + stream1.str() + "cm , |ydiff| < " + stream2.str() + "cm"),
0198                                       std::string("xy cuts and |etadiff| < " + stream3.str()),
0199                                       std::string("xy, eta cuts and |phidiff| < " + stream4.str()),
0200                                       "All Track Seeds (North Only)",
0201                                       "All Track Seeds (South Only)",
0202                                       "North All Cuts (x,y,eta,phi)",
0203                                       "South All Cuts (x,y,eta,phi)"};
0204
0205   {
0206     h_crossing = new TH1F(std::string(getHistoPrefix() + "crossing").c_str(),
0207                           "Track Crossing Value", 1000, -500, 500);
0208     h_crossing->GetXaxis()->SetTitle("Track Crossing");
0209     h_crossing->GetYaxis()->SetTitle("Entries");
0210     hm->registerHisto(h_crossing);
0211   }
0212   int i = 0;
0213   for (const std::string& name : cutNames)
0214   {
0215     std::string histoTitle = std::string(cutVals[i]);
0216     h_phiDiff[i] = new TH1F(std::string(getHistoPrefix() + "phiDiff" + name).c_str(),
0217                             histoTitle.c_str(), 100, -0.5, 0.5);
0218     h_phiDiff[i]->GetXaxis()->SetTitle("TPC Seed #phi - Silicon Seed #phi");
0219     h_phiDiff[i]->GetYaxis()->SetTitle("Entries");
0220     hm->registerHisto(h_phiDiff[i]);
0221     i++;
0222   }
0223   i = 0;
0224   for (const std::string& name : cutNames)
0225   {
0226     std::string histoTitle = std::string(cutVals[i]);
0227     h_etaDiff[i] = new TH1F(std::string(getHistoPrefix() + "etaDiff" + name).c_str(),
0228                             histoTitle.c_str(), 100, -0.1, 0.1);
0229     h_etaDiff[i]->GetXaxis()->SetTitle("TPC Seed #eta - Silicon Seed #eta");
0230     h_etaDiff[i]->GetYaxis()->SetTitle("Entries");
0231     hm->registerHisto(h_etaDiff[i]);
0232     i++;
0233   }
0234   i = 0;
0235   for (const std::string& name : cutNames)
0236   {
0237     std::string histoTitle = std::string(cutVals[i]);
0238     h_xDiff[i] = new TH1F(std::string(getHistoPrefix() + "xDiff" + name).c_str(),
0239                           histoTitle.c_str(), 100, -2, 2);
0240     h_xDiff[i]->GetXaxis()->SetTitle("TPC Seed x - Silicon Seed x [cm]");
0241     h_xDiff[i]->GetYaxis()->SetTitle("Entries");
0242     hm->registerHisto(h_xDiff[i]);
0243     i++;
0244   }
0245   i = 0;
0246   for (const std::string& name : cutNames)
0247   {
0248     std::string histoTitle = std::string(cutVals[i]);
0249     h_yDiff[i] = new TH1F(std::string(getHistoPrefix() + "yDiff" + name).c_str(),
0250                           histoTitle.c_str(), 100, -2, 2);
0251     h_yDiff[i]->GetXaxis()->SetTitle("TPC Seed y - Silicon Seed y [cm]");
0252     h_yDiff[i]->GetYaxis()->SetTitle("Entries");
0253     hm->registerHisto(h_yDiff[i]);
0254     i++;
0255   }
0256   i = 0;
0257   for (const std::string& name : cutNames)
0258   {
0259     std::string histoTitle = std::string(cutVals[i]);
0260     h_zDiff[i] = new TH1F(std::string(getHistoPrefix() + "zDiff" + name).c_str(),
0261                           histoTitle.c_str(), 500, -100, 100);
0262     h_zDiff[i]->GetXaxis()->SetTitle("TPC Seed z - Silicon Seed z [cm]");
0263     h_zDiff[i]->GetYaxis()->SetTitle("Entries");
0264     hm->registerHisto(h_zDiff[i]);
0265     i++;
0266   }
0267
0268   return;
0269 }