tpc/tpcchannelmasking/TpcMaskedChannelMap.cc

0001
0002 #include "TpcMaskedChannelMap.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <fun4all/Fun4AllServer.h>
0006
0007 #include <phool/getClass.h>
0008
0009 #include <ffarawobjects/TpcRawHit.h>
0010 #include <ffarawobjects/TpcRawHitContainer.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>    // for PHIODataNode
0014 #include <phool/PHNode.h>          // for PHNode
0015 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0016 #include <phool/PHPointerListIterator.h>
0017
0018 #include <cdbobjects/CDBTTree.h>
0019 #include <g4detectors/PHG4TpcGeom.h>
0020 #include <g4detectors/PHG4TpcGeomContainer.h>
0021 #include <trackbase/TpcDefs.h>
0022
0023 #include <TFile.h>
0024
0025 #include <cassert>
0026 #include <iostream>
0027 #include <memory>
0028 #include <regex>
0029
0030
0031 //____________________________________________________________________________..
0032 TpcMaskedChannelMap::TpcMaskedChannelMap(const std::string &name)
0033   : SubsysReco("TpcMaskedChannelMap")
0034   , m_run(name)
0035 {
0036   M.setMapNames("AutoPad-R1-RevA.sch.ChannelMapping.csv", "AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv", "AutoPad-R3-RevA.sch.ChannelMapping.csv");
0037 }
0038
0039 //____________________________________________________________________________..
0040 int TpcMaskedChannelMap::InitRun(PHCompositeNode *topNode)
0041 {
0042   PHNodeIterator trkr_itr(topNode);
0043   PHCompositeNode *tpc_node = dynamic_cast<PHCompositeNode *>(trkr_itr.findFirst("PHCompositeNode", "TPC"));
0044   if (!tpc_node)
0045   {
0046     std::cout << __PRETTY_FUNCTION__ << " : ERROR : "
0047               << "No TPC node found, unregistering module" << std::endl;
0048     Fun4AllServer::instance()->unregisterSubsystem(this);
0049     return Fun4AllReturnCodes::EVENT_OK;
0050   }
0051
0052   PHNodeIterator tpc_itr(tpc_node);
0053   {
0054     PHPointerListIterator<PHNode> iter(tpc_itr.ls());
0055
0056     PHNode *thisNode_raw;
0057     while ((thisNode_raw = iter()))
0058     {
0059       if (thisNode_raw->getType() != "PHIODataNode")
0060       {
0061         continue;
0062       }
0063
0064       PHIODataNode<TpcRawHitContainer> *thisNode = static_cast<PHIODataNode<TpcRawHitContainer> *>(thisNode_raw);
0065       if (thisNode)
0066       {
0067         std::cout << __PRETTY_FUNCTION__ << " : Found TpcRawHitContainer Node "
0068                   << thisNode->getName() << std::endl;
0069
0070         TpcRawHitContainer *rawhitcont = (TpcRawHitContainer *) thisNode->getData();
0071         if (rawhitcont)
0072         {
0073           rawhitcont_vec.push_back(rawhitcont);
0074         }
0075       }
0076     }
0077   }
0078
0079   m_hotFile = "TPCHotMap_" + m_run + ".root";
0080   m_deadFile = "TPCDeadMap_" + m_run + ".root";
0081
0082   m_histogramFile = "HIST_" + m_run + ".root";
0083
0084   h_hits_side0 = new TH2F("h_hits_side0","Hits (Side0);X [cm];Y [cm]",400,-800,800,400,-800,800);
0085   h_hits_side1 = new TH2F("h_hits_side1","Hits (Side1);X [cm];Y [cm]",400,-800,800,400,-800,800);
0086   h_masked_side0 = new TH2F("h_masked_side0","Masked Channels (Side0);X [cm];Y [cm]",400,-800,800,400,-800,800);
0087   h_masked_side1 = new TH2F("h_masked_side1","Masked Channels (Side1);X [cm];Y [cm]",400,-800,800,400,-800,800);
0088
0089   return Fun4AllReturnCodes::EVENT_OK;
0090 }
0091
0092 //____________________________________________________________________________..
0093 int TpcMaskedChannelMap::process_event(PHCompositeNode * /* unused */)
0094 {
0095   n_Events += 1.;
0096
0097   unsigned int raw_hit_num = 0;
0098   for (TpcRawHitContainer *&rawhitcont : rawhitcont_vec)
0099   {
0100     raw_hit_num = rawhitcont->get_nhits();
0101     for (unsigned int i = 0; i < raw_hit_num; i++)
0102     {
0103       auto *hit = rawhitcont->get_hit(i);
0104       int32_t packet_id = hit->get_packetid();
0105       int ep = (packet_id - 4000) % 10;
0106       int sector = (packet_id - 4000 - ep) / 10;
0107       uint16_t fee = hit->get_fee();
0108       int channel = hit->get_channel();
0109
0110       if (fee == 16 && channel < 32)
0111       {
0112         continue;
0113       }
0114
0115       int feeM = FEE_map[fee];
0116       if (FEE_R[fee] == 2)
0117       {
0118         feeM += 6;
0119       }
0120       if (FEE_R[fee] == 3)
0121       {
0122         feeM += 14;
0123       }
0124       double R = M.getR(feeM, channel);
0125       double phi = 0;
0126       if (sector < 12)  // NS
0127       {
0128         phi = M.getPhi(feeM, channel) + (sector) *M_PI / 6;
0129       }
0130       else if (sector >= 12)  // SS
0131       {
0132         phi = M.getPhi(feeM, channel) + (sector-12) * M_PI / 6;
0133       }
0134
0135       float median = 60;
0136       for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(hit->CreateAdcIterator());
0137            !adc_iterator->IsDone();
0138            adc_iterator->Next())
0139       {
0140         const uint16_t sampleN = adc_iterator->CurrentTimeBin();
0141         const uint16_t adc = adc_iterator->CurrentAdc();
0142         if (adc - median <= 20)
0143         {
0144           continue;
0145         }
0146
0147         if (sampleN >= 400 && sampleN <= 430)
0148         {
0149           continue;
0150         }
0151
0152         nhit_sectors_fees_channels[sector][feeM][channel] += 1;
0153
0154         if (sector < 12)
0155         {
0156           h_hits_side1->Fill(R * cos(phi), R * sin(phi));
0157         }
0158         else
0159         {
0160           h_hits_side0->Fill(R * cos(phi), R * sin(phi));
0161         }
0162       }
0163     }
0164   }
0165
0166   if (raw_hit_num == 0)
0167   {
0168     return Fun4AllReturnCodes::EVENT_OK;
0169   }
0170
0171   return Fun4AllReturnCodes::EVENT_OK;
0172 }
0173
0174 //____________________________________________________________________________..
0175 int TpcMaskedChannelMap::End(PHCompositeNode *topNode)
0176 {
0177   PHG4TpcGeomContainer* geom_container = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0178   if (!geom_container)
0179   {
0180     std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << std::endl;
0181     return Fun4AllReturnCodes::ABORTRUN;
0182   }
0183
0184   int nDeadFees = 0;
0185
0186   int side = 1;
0187   for (int m_Sector = 0; m_Sector < 24; m_Sector++)
0188   {
0189     if (m_Sector > 11)
0190     {
0191       side = 0;
0192     }
0193     for (int m_FEE = 0; m_FEE < 26; m_FEE++)
0194     {
0195       int nDeadInFee = 0;
0196       for (int m_Channel = 0; m_Channel < 256; m_Channel++)
0197       {
0198         if (nhit_sectors_fees_channels[m_Sector][m_FEE][m_Channel]/n_Events < m_deadChanHitCut)
0199         {
0200           nDeadInFee++;
0201           int layer = M.getLayer(m_FEE, m_Channel);
0202
0203           if (layer < 7)
0204           {
0205             continue;
0206           }
0207
0208           PHG4TpcGeom* layergeom = geom_container->GetLayerCellGeom(layer);
0209
0210           double phi = 0;
0211           if (m_Sector < 12)  // NS
0212           {
0213             phi = M.getPhi(m_FEE, m_Channel) + (m_Sector) *M_PI / 6;
0214           }
0215           else if (m_Sector >= 12)  // SS
0216           {
0217             phi = M.getPhi(m_FEE, m_Channel) + (m_Sector-12) * M_PI / 6;
0218           }
0219           unsigned int phibin = layergeom->find_phibin(phi,side);
0220           double R = M.getR(m_FEE, m_Channel);
0221           m_deadChannelCDB.emplace(layer, mc_sectors[m_Sector % 12], side, phibin,R * cos(phi), R * sin(phi));
0222           if (side == 0)
0223           {
0224             h_masked_side0->Fill(R * cos(phi), R * sin(phi));
0225           }
0226           else
0227           {
0228             h_masked_side1->Fill(R * cos(phi), R * sin(phi));
0229           }
0230         }
0231         else if (nhit_sectors_fees_channels[m_Sector][m_FEE][m_Channel]/n_Events > m_hotChanHitCut)
0232         {
0233           nDeadInFee++;
0234           int layer = M.getLayer(m_FEE, m_Channel);
0235
0236           if (layer < 7)
0237           {
0238             continue;
0239           }
0240
0241           PHG4TpcGeom* layergeom = geom_container->GetLayerCellGeom(layer);
0242
0243           double phi = 0;
0244           if (m_Sector < 12)  // NS
0245           {
0246             phi = M.getPhi(m_FEE, m_Channel) + (m_Sector) *M_PI / 6;
0247           }
0248           else if (m_Sector >= 12)  // SS
0249           {
0250             phi = M.getPhi(m_FEE, m_Channel) + (m_Sector-12) * M_PI / 6;
0251           }
0252           unsigned int phibin = layergeom->find_phibin(phi,side);
0253           double R = M.getR(m_FEE, m_Channel);
0254           m_hotChannelCDB.emplace(layer, mc_sectors[m_Sector % 12], side, phibin,R * cos(phi), R * sin(phi));
0255           if (side == 0)
0256           {
0257             h_masked_side0->Fill(R * cos(phi), R * sin(phi));
0258           }
0259           else
0260           {
0261             h_masked_side1->Fill(R * cos(phi), R * sin(phi));
0262           }
0263         }
0264       }
0265       if (nDeadInFee > 250)
0266       {
0267         nDeadFees++;
0268       }
0269     }
0270   }
0271
0272   std::cout << "Number of dead FEEs: " << nDeadFees << std::endl;
0273
0274   CDBTTree cdbttree( m_deadFile );
0275   cdbttree.SetSingleIntValue( "TotalDeadChannels", m_deadChannelCDB.size() );
0276
0277   std::cout << "Total Number of Dead Channels: " << m_deadChannelCDB.size() << std::endl;
0278
0279   int index = 0;
0280   for( const auto& channel_id:m_deadChannelCDB )
0281   {
0282     cdbttree.SetIntValue( index, "layer", channel_id.m_layer );
0283     cdbttree.SetIntValue( index, "sector", channel_id.m_sector );
0284     cdbttree.SetIntValue( index, "side", channel_id.m_side );
0285     cdbttree.SetIntValue( index, "pad", channel_id.m_pad );
0286     cdbttree.SetFloatValue( index, "x", channel_id.m_x );
0287     cdbttree.SetFloatValue( index, "y", channel_id.m_y );
0288     ++index;
0289   }
0290
0291   // commit and write
0292   cdbttree.Commit();
0293   cdbttree.CommitSingle();
0294   cdbttree.WriteCDBTTree();
0295   std::cout << __PRETTY_FUNCTION__ << " : completed saving to " << m_deadFile << std::endl;
0296
0297   CDBTTree cdbttree2( m_hotFile );
0298   cdbttree2.SetSingleIntValue( "TotalHotChannels", m_hotChannelCDB.size() );
0299
0300   std::cout << "Total Number of Hot Channels: " << m_hotChannelCDB.size() << std::endl;
0301
0302   index = 0;
0303   for( const auto& channel_id:m_hotChannelCDB )
0304   {
0305     cdbttree2.SetIntValue( index, "layer", channel_id.m_layer );
0306     cdbttree2.SetIntValue( index, "sector", channel_id.m_sector );
0307     cdbttree2.SetIntValue( index, "side", channel_id.m_side );
0308     cdbttree2.SetIntValue( index, "pad", channel_id.m_pad );
0309     cdbttree2.SetFloatValue( index, "x", channel_id.m_x );
0310     cdbttree2.SetFloatValue( index, "y", channel_id.m_y );
0311     ++index;
0312   }
0313
0314   // commit and write
0315   cdbttree2.Commit();
0316   cdbttree2.CommitSingle();
0317   cdbttree2.WriteCDBTTree();
0318   std::cout << __PRETTY_FUNCTION__ << " : completed saving to " << m_hotFile << std::endl;
0319
0320   TFile *outFile = new TFile(m_histogramFile.c_str(), "RECREATE");
0321   h_hits_side0->Write();
0322   h_hits_side1->Write();
0323   h_masked_side0->Write();
0324   h_masked_side1->Write();
0325   outFile->Close();
0326
0327   delete h_hits_side0;
0328   delete h_hits_side1;
0329   delete h_masked_side0;
0330   delete h_masked_side1;
0331   delete outFile;
0332
0333   return Fun4AllReturnCodes::EVENT_OK;
0334 }