Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:20:32 

0001 #include "TpcRawWriter.h"
0002 
0003 #include <micromegas/MicromegasDefs.h>
0004 #include <trackbase/InttDefs.h>
0005 #include <trackbase/MvtxDefs.h>
0006 #include <trackbase/TpcDefs.h>
0007 
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009 #include <fun4all/SubsysReco.h>  // for SubsysReco
0010 #include <trackbase/ActsGeometry.h>
0011 #include <trackbase/RawHitSet.h>
0012 #include <trackbase/RawHitSetContainer.h>
0013 #include <trackbase/RawHitSetContainerv1.h>
0014 #include <trackbase/RawHitSetv1.h>
0015 #include <trackbase/RawHitTpc.h>
0016 #include <trackbase/RawHitv1.h>
0017 #include <trackbase/TrkrCluster.h>
0018 #include <trackbase/TrkrClusterContainerv4.h>
0019 #include <trackbase/TrkrClusterHitAssocv3.h>
0020 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0021 #include <trackbase/TrkrHit.h>
0022 #include <trackbase/TrkrHitSet.h>
0023 #include <trackbase/TrkrHitSetContainer.h>
0024 
0025 #include <g4detectors/PHG4TpcGeom.h>
0026 #include <g4detectors/PHG4TpcGeomContainer.h>
0027 
0028 #include <Acts/Definitions/Units.hpp>
0029 #include <Acts/Surfaces/Surface.hpp>
0030 
0031 #include <phool/PHCompositeNode.h>
0032 #include <phool/PHIODataNode.h>  // for PHIODataNode
0033 #include <phool/PHNode.h>        // for PHNode
0034 #include <phool/PHNodeIterator.h>
0035 #include <phool/PHObject.h>  // for PHObject
0036 #include <phool/getClass.h>
0037 #include <phool/phool.h>  // for PHWHERE
0038 
0039 #include <TMatrixFfwd.h>    // for TMatrixF
0040 #include <TMatrixT.h>       // for TMatrixT, ope...
0041 #include <TMatrixTUtils.h>  // for TMatrixTRow
0042 
0043 #include <TFile.h>
0044 
0045 #include <array>
0046 #include <cmath>  // for sqrt, cos, sin
0047 #include <iostream>
0048 #include <limits>
0049 #include <map>  // for _Rb_tree_cons...
0050 #include <string>
0051 #include <utility>  // for pair
0052 #include <vector>
0053 // Terra incognita....
0054 #include <pthread.h>
0055 
0056 TpcRawWriter::TpcRawWriter(const std::string &name)
0057   : SubsysReco(name)
0058 {
0059   std::cout << PHWHERE << "Construct TpcRawWriter" << std::endl;
0060 }
0061 
0062 int TpcRawWriter::InitRun(PHCompositeNode *topNode)
0063 {
0064   if (topNode)
0065   {
0066     std::cout << PHWHERE << "Init TpcRawWriter" << std::endl;
0067   }
0068   PHNodeIterator iter(topNode);
0069 
0070   // Looking for the DST node
0071   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0072   if (!dstNode)
0073   {
0074     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0075     return Fun4AllReturnCodes::ABORTRUN;
0076   }
0077 
0078   // Create the Cluster node if required
0079   auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
0080   if (!trkrclusters)
0081   {
0082     PHNodeIterator dstiter(dstNode);
0083     PHCompositeNode *DetNode =
0084         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0085     if (!DetNode)
0086     {
0087       DetNode = new PHCompositeNode("TRKR");
0088       dstNode->addNode(DetNode);
0089     }
0090 
0091     trkrclusters = new TrkrClusterContainerv4;
0092     PHIODataNode<PHObject> *TrkrClusterContainerNode =
0093         new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
0094     DetNode->addNode(TrkrClusterContainerNode);
0095   }
0096 
0097   auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0098   if (!clusterhitassoc)
0099   {
0100     PHNodeIterator dstiter(dstNode);
0101     PHCompositeNode *DetNode =
0102         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0103     if (!DetNode)
0104     {
0105       DetNode = new PHCompositeNode("TRKR");
0106       dstNode->addNode(DetNode);
0107     }
0108 
0109     clusterhitassoc = new TrkrClusterHitAssocv3;
0110     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
0111     DetNode->addNode(newNode);
0112   }
0113 
0114   // new containers
0115   m_rawhits = findNode::getClass<RawHitSetContainerv1>(topNode, "TRKR_RAWHITSET");
0116   if (!m_rawhits)
0117   {
0118     PHNodeIterator dstiter(dstNode);
0119     auto DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0120     if (!DetNode)
0121     {
0122       DetNode = new PHCompositeNode("TRKR");
0123       dstNode->addNode(DetNode);
0124     }
0125 
0126     m_rawhits = new RawHitSetContainerv1;
0127     auto newNode = new PHIODataNode<PHObject>(m_rawhits, "TRKR_RAWHITSET", "PHObject");
0128     DetNode->addNode(newNode);
0129   }
0130 
0131   return Fun4AllReturnCodes::EVENT_OK;
0132 }
0133 
0134 int TpcRawWriter::process_event(PHCompositeNode *topNode)
0135 {
0136   //  int print_layer = 18;
0137 
0138   //  if (Verbosity() > 1000)
0139   if (topNode)
0140   {
0141     std::cout << "TpcRawWriter::Process_Event" << std::endl;
0142   }
0143 
0144   PHNodeIterator iter(topNode);
0145   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0146   if (!dstNode)
0147   {
0148     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0149     return Fun4AllReturnCodes::ABORTRUN;
0150   }
0151 
0152   // get node containing the digitized hits
0153   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0154   if (!m_hits)
0155   {
0156     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0157     return Fun4AllReturnCodes::ABORTRUN;
0158   }
0159 
0160   // new containers
0161   m_rawhits = findNode::getClass<RawHitSetContainerv1>(topNode, "TRKR_RAWHITSET");
0162   if (!m_rawhits)
0163   {
0164     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0165     return Fun4AllReturnCodes::ABORTRUN;
0166   }
0167 
0168   // get node for clusters
0169   m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0170   if (!m_clusterlist)
0171   {
0172     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
0173     return Fun4AllReturnCodes::ABORTRUN;
0174   }
0175 
0176   // get node for cluster hit associations
0177   m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0178   if (!m_clusterhitassoc)
0179   {
0180     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
0181     return Fun4AllReturnCodes::ABORTRUN;
0182   }
0183 
0184   PHG4TpcGeomContainer *geom_container =
0185       findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0186   if (!geom_container)
0187   {
0188     std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << std::endl;
0189     return Fun4AllReturnCodes::ABORTRUN;
0190   }
0191 
0192   m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
0193                                                  "ActsGeometry");
0194   if (!m_tGeometry)
0195   {
0196     std::cout << PHWHERE
0197               << "ActsTrackingGeometry not found on node tree. Exiting"
0198               << std::endl;
0199     return Fun4AllReturnCodes::ABORTRUN;
0200   }
0201 
0202   // The hits are stored in hitsets, where each hitset contains all hits in a given TPC readout (layer, sector, side), so clusters are confined to a hitset
0203   // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
0204 
0205   // loop over the MVTX HitSet objects
0206   std::cout << "processing mvtx" << std::endl;
0207   TrkrHitSetContainer::ConstRange mvtxhitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::mvtxId);
0208   //  const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
0209 
0210 //  int count = 0;
0211   for (TrkrHitSetContainer::ConstIterator hitsetitr = mvtxhitsetrange.first;
0212        hitsetitr != mvtxhitsetrange.second;
0213        ++hitsetitr)
0214   {
0215     TrkrHitSet *hitset = hitsetitr->second;
0216 
0217     m_rawhits->findOrAddHitSet(hitsetitr->first);
0218     RawHitSetv1 *rhitset = dynamic_cast<RawHitSetv1 *>(m_rawhits->findHitSet(hitsetitr->first));
0219     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0220 
0221     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0222          hitr != hitrangei.second;
0223          ++hitr)
0224     {
0225       unsigned short adc = (unsigned short) (hitr->second->getAdc());
0226       if (adc > 0)
0227       {
0228         RawHitv1 *rHit = new RawHitv1;
0229         unsigned short iphi = MvtxDefs::getCol(hitr->first);
0230         unsigned short it = MvtxDefs::getRow(hitr->first);
0231         rHit->setPhiBin(iphi);
0232         rHit->setTBin(it);
0233         rHit->setAdc(adc);
0234         rhitset->addHit(rHit);
0235       }
0236     }
0237 
0238 //    count++;
0239   }
0240   std::cout << "processing intt" << std::endl;
0241   // loop over the INTT HitSet objects
0242   TrkrHitSetContainer::ConstRange intt_hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::inttId);
0243   //  const int num_hitsets = std::distance(intt_hitsetrange.first,intt_hitsetrange.second);
0244 
0245   for (TrkrHitSetContainer::ConstIterator hitsetitr = intt_hitsetrange.first;
0246        hitsetitr != intt_hitsetrange.second;
0247        ++hitsetitr)
0248   {
0249     TrkrHitSet *hitset = hitsetitr->second;
0250 
0251     m_rawhits->findOrAddHitSet(hitsetitr->first);
0252     RawHitSet *rhitset = m_rawhits->findHitSet(hitsetitr->first);
0253     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0254 
0255     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0256          hitr != hitrangei.second;
0257          ++hitr)
0258     {
0259       // int sector = InttDefs::getLadderPhiId(hitsetitr->first);
0260       // int side = InttDefs::getLadderZId(hitsetitr->first);
0261       // int layer  = TrkrDefs::getLayer(hitsetitr->first);
0262       unsigned short adc = (unsigned short) (hitr->second->getAdc());
0263       if (adc > 0)
0264       {
0265         RawHitv1 *rHit = new RawHitv1;
0266         unsigned short iphi = InttDefs::getCol(hitr->first);
0267         unsigned short it = InttDefs::getRow(hitr->first);
0268         // std::cout << " layer " << layer << " sector: " << sector << " side " << side << " col: " << iphi << " row " << it << std::endl;
0269         rHit->setPhiBin(iphi);
0270         rHit->setTBin(it);
0271         rHit->setAdc(adc);
0272         rhitset->addHit(rHit);
0273       }
0274     }
0275 
0276 //    count++;
0277   }
0278   std::cout << "processing mm" << std::endl;
0279   // loop over the micromega HitSet objects
0280   TrkrHitSetContainer::ConstRange mm_hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::micromegasId);
0281   //  const int num_hitsets = std::distance(mm_hitsetrange.first,mm_hitsetrange.second);
0282 
0283   for (TrkrHitSetContainer::ConstIterator hitsetitr = mm_hitsetrange.first;
0284        hitsetitr != mm_hitsetrange.second;
0285        ++hitsetitr)
0286   {
0287     TrkrHitSet *hitset = hitsetitr->second;
0288 
0289     m_rawhits->findOrAddHitSet(hitsetitr->first);
0290     RawHitSet *rhitset = m_rawhits->findHitSet(hitsetitr->first);
0291     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0292 
0293     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0294          hitr != hitrangei.second;
0295          ++hitr)
0296     {
0297       unsigned short adc = (unsigned short) (hitr->second->getAdc());
0298       if (adc > 0)
0299       {
0300         RawHitv1 *rHit = new RawHitv1;
0301         //  TrkrDefs::hitkey tmp = (hitr->first >> MicromegasDefs::kBitShiftStrip);
0302         unsigned short iphi = MicromegasDefs::getStrip(hitr->first);
0303         unsigned short it = 0;
0304         rHit->setPhiBin(iphi);
0305         rHit->setTBin(it);
0306         rHit->setAdc(adc);
0307         rhitset->addHit(rHit);
0308       }
0309     }
0310 
0311 //    count++;
0312   }
0313   std::cout << "processing tpc" << std::endl;
0314   float tpc_zmax = m_tGeometry->get_max_driftlength() + m_tGeometry->get_CM_halfwidth();
0315 
0316   // loop over the TPC HitSet objects
0317   TrkrHitSetContainer::ConstRange tpc_hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0318   //  const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
0319 
0320   int ncheck = 0;
0321   int allbad = 0;
0322 
0323   for (TrkrHitSetContainer::ConstIterator hitsetitr = tpc_hitsetrange.first;
0324        hitsetitr != tpc_hitsetrange.second;
0325        ++hitsetitr)
0326   // TrkrHitSetContainer::ConstIterator hitsetitr = tpc_hitsetrange.first;
0327   {
0328     TrkrHitSet *hitset = hitsetitr->second;
0329     unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0330     //   if(layer!=7) continue;
0331     //    int side = TpcDefs::getSide(hitsetitr->first);
0332     unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0333     PHG4TpcGeom *layergeom = geom_container->GetLayerCellGeom(layer);
0334     unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
0335     unsigned short NPhiBinsSector = NPhiBins / 12;
0336     unsigned short NZBins = (unsigned short) layergeom->get_zbins();
0337     // unsigned short NZBinsSide = NZBins/2;
0338     unsigned short NZBinsSide = NZBins;
0339     unsigned short NZBinsMin = 0;
0340     unsigned short PhiOffset = NPhiBinsSector * sector;
0341     unsigned short ZOffset = NZBinsMin;
0342     std::vector<int> nhits_thispad;
0343     std::vector<std::vector<unsigned short>> adcval;
0344     adcval.resize(NPhiBins);
0345     for (int nphi = 0; nphi < NPhiBins; nphi++)
0346     {
0347       for (int nz = 0; nz < NZBins; nz++)
0348       {
0349         adcval[nphi].push_back(0);
0350         nhits_thispad.push_back(0);
0351       }
0352     }
0353     // std::vector<std::vector<unsigned short>> outadc;
0354     // outadc.resize(NPhiBins);
0355 
0356     //    RawHitSetContainerv1::Iterator rawhitset_iter =
0357     m_rawhits->findOrAddHitSet(hitsetitr->first);
0358     RawHitSetv1 *rhitset = dynamic_cast<RawHitSetv1 *>(m_rawhits->findHitSet(hitsetitr->first));
0359     rhitset->setTpcPhiBins(NPhiBins);
0360 
0361     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0362     int zbinmax = 498;
0363     int zbinmin = 0;
0364     if (layer >= 7 && layer < 22)
0365     {
0366       int etacut = 249 - ((50 + (layer - 7)) / tpc_zmax) * 249;
0367       zbinmin = etacut;
0368       zbinmax -= etacut;
0369     }
0370     if (layer >= 22 && layer <= 48)
0371     {
0372       int etacut = 249 - ((65 + ((40.5 / 26) * (layer - 22))) / tpc_zmax) * 249;
0373       zbinmin = etacut;
0374       zbinmax -= etacut;
0375     }
0376 
0377     //    std::cout << " layer: " << layer << " zbin limit " << zbinmin << " | " << zbinmax <<std::endl;
0378     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0379          hitr != hitrangei.second;
0380          ++hitr)
0381     {
0382       if (TpcDefs::getPad(hitr->first) - PhiOffset < 0)
0383       {
0384         continue;
0385       }
0386       if (TpcDefs::getTBin(hitr->first) - ZOffset < 0)
0387       {
0388         // std::cout << "WARNING zbin out of range: " << TpcDefs::getTBin(hitr->first) - zoffset  << " | " << zbins <<std::endl;
0389         continue;
0390       }
0391       unsigned short phibin = TpcDefs::getPad(hitr->first) - PhiOffset;
0392       unsigned short zbin = TpcDefs::getTBin(hitr->first) - ZOffset;
0393       unsigned short zbinorg = TpcDefs::getTBin(hitr->first);
0394       if (phibin >= NPhiBinsSector)
0395       {
0396         // std::cout << "WARNING phibin out of range: " << phibin << " | " << phibins << std::endl;
0397         continue;
0398       }
0399       if (zbin >= NZBinsSide)
0400       {
0401         // std::cout << "WARNING z bin out of range: " << zbin << " | " << zbins << std::endl;
0402         continue;
0403       }
0404       if (zbinorg > zbinmax || zbinorg < zbinmin)
0405       {
0406         continue;
0407       }
0408       float_t fadc = (hitr->second->getAdc()) - pedestal;  // proper int rounding +0.5
0409       unsigned short adc = 0;
0410       if (fadc > 0)
0411       {
0412         adc = (unsigned short) fadc;
0413       }
0414 
0415       if (adc > 0)
0416       {
0417         nhits_thispad[phibin]++;
0418         unsigned short it = TpcDefs::getTBin(hitr->first);
0419         adcval[phibin][it] = adc;
0420         // if(layer==20&&sector==4&&phibin==77)
0421         // std::cout << "input sector: " << sector << " lay: " << layer << " # " << phibin << "|"<< it <<" adc: " << (int)adcval[phibin][it] << " pack: " << std::endl;
0422         /*RawHitTpc  *rHit = new RawHitTpc;
0423         unsigned short it = TpcDefs::getTBin(hitr->first);
0424         rHit->setTBin(it);
0425         unsigned short iadc = adc;
0426         rHit->setAdc(iadc);
0427         rhitset->addTpcHit(phibin,rHit);
0428         */
0429       }
0430     }
0431     //    assert(adcval);
0432     for (int nphi = 0; nphi < NPhiBins; nphi++)
0433     {
0434       for (int nz = 0; nz < NZBins; nz++)
0435       {
0436         /*if(nphi>=1&&nz<NZBins-1)//Remove isolated single timebin hits
0437           if(layer==29&&sector==0&&nphi==0){
0438             if(nz>209&&nz<215){
0439               std::cout << "single check sector: " << sector << " lay: " << layer << " # " << nphi << "|"<< nz <<" adc: " << (int)adcval[nphi][nz] << " pack: " << std::endl;
0440             }
0441           }
0442         */
0443         if (adcval[nphi][nz] > 0)
0444         {
0445           if ((adcval[nphi][nz - 1] == 0) && (adcval[nphi][nz + 1] == 0))
0446           {
0447             // std::cout << "single hit sector: " << sector << " lay: " << layer << " # " << nphi << "|"<< nz <<" adc: " << (int)adcval[nphi][nz] << " pack: " << std::endl;
0448             adcval[nphi][nz] = 0;
0449           }
0450         }
0451       }
0452     }
0453     for (int nphi = 0; nphi < NPhiBins; nphi++)
0454     {
0455       int zero_count = 0;
0456       // int outpos = 0;
0457       if (nhits_thispad[nphi] == 0)
0458       {
0459         continue;
0460       }
0461       for (int nz = 0; nz < NZBins; nz++)
0462       {
0463         uint8_t thisadc = 0;
0464         if (nphi >= 1 && nz < NZBins - 1)
0465         {  // Remove isolated single timebin hits
0466           if (adcval[nphi][nz] > 0)
0467           {
0468             if ((adcval[nphi][nz - 1] == 0) && (adcval[nphi][nz + 1] == 0))
0469             {
0470               adcval[nphi][nz] = 0;
0471             }
0472           }
0473         }
0474 
0475         if (adcval[nphi][nz] > 255)
0476         {  // take care of overflows, we want to store 8 bit ADC values
0477           thisadc = 255;
0478           //        adcval[nphi][nz]=255;
0479         }
0480         else
0481         {
0482           thisadc = adcval[nphi][nz];
0483         }
0484         // std::cout << "nz: " << nz << " adc " << (int)thisadc << std::endl;
0485         if (thisadc > 0)
0486         {  // non zero adc, fill adc
0487           zero_count = 0;
0488           (*(rhitset->getHits(nphi))).push_back(thisadc);
0489           if (layer == 222 && sector == 6 && nphi == 10)
0490           {
0491             std::cout << "0#" << nphi << "nz= " << nz << " filling " << (int) thisadc << " at " << rhitset->size(nphi) << std::endl;
0492           }
0493         }
0494         else
0495         {
0496           zero_count++;
0497           //      if(nz==0){
0498           //  rhitset->m_tpchits[nphi].push_back(0);
0499           // std::cout << " 1filling " << 0 << "|" << rhitset->m_tpchits[nphi].back() << " | " << rhitset->m_tpchits[nphi][outpos++] << std::endl;
0500           // if(layer==222&&sector==6&&nphi==10)std::cout << "1#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
0501           // }else
0502           /*
0503           if(adcval[nphi][nz-1]>0){ //trailing edge 0
0504             rhitset->m_tpchits[nphi].push_back(0);
0505             if(layer==222&&sector==6&&nphi==10)std::cout << "2#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
0506             //std::cout << " 2filling " << 0 << "|" << rhitset->m_tpchits[nphi].back() << " | " << rhitset->m_tpchits[nphi][outpos++] << std::endl;
0507           }
0508           */
0509         }
0510         if (zero_count > 1)
0511         {
0512           if (nz < NZBins - 2)
0513           {
0514             if (adcval[nphi][nz + 2] > 0 && adcval[nphi][nz + 1] == 0)
0515             {  // fill zero count, end of zero series
0516               (*(rhitset->getHits(nphi))).push_back(zero_count - 1);
0517               zero_count = 0;
0518               if (layer == 222 && sector == 6 && nphi == 10)
0519               {
0520                 std::cout << "3#" << nphi << "nz= " << nz << " filling " << zero_count - 1 << " at " << rhitset->size(nphi) << " nz+2 " << (nz + 2) << " adcval[nz+2]" << (int) adcval[nphi][nz + 2] << std::endl;
0521               }
0522             }
0523             if (adcval[nphi][nz + 1] > 0)
0524             {
0525               (*(rhitset->getHits(nphi))).push_back(0);  // leading edge zero
0526               if (layer == 222 && sector == 6 && nphi == 10)
0527               {
0528                 std::cout << "4#" << nphi << "nz= " << nz << " filling " << 0 << " at " << rhitset->size(nphi) << std::endl;
0529               }
0530 
0531               // std::cout << " 4filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
0532             }
0533           }
0534         }
0535         else
0536         {
0537           if (zero_count == 1)
0538           {
0539             (*(rhitset->getHits(nphi))).push_back(0);
0540             if (layer == 222 && sector == 6 && nphi == 10)
0541             {
0542               std::cout << "5#" << nphi << "nz= " << nz << " filling " << 0 << " at " << rhitset->size(nphi) << std::endl;
0543             }
0544             // std::cout << " 5filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
0545           }
0546         }
0547         if (zero_count == 254)
0548         {
0549           zero_count = 0;
0550           (*(rhitset->getHits(nphi))).push_back(254 - 2);
0551           (*(rhitset->getHits(nphi))).push_back(0);
0552           if (layer == 222 && sector == 6 && nphi == 10)
0553           {
0554             std::cout << "6#" << nphi << "nz= " << nz << " filling " << 254 - 1 << " at " << rhitset->size(nphi) << std::endl;
0555           }
0556           // std::cout << "6 filling " << 254-1 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
0557           //  rhitset->m_tpchits[nphi].push_back(0);
0558           //  if(layer==222&&sector==6&&nphi==10)std::cout << "7#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
0559           // std::cout << " 7filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
0560         }
0561         /*
0562           if(layer%15==0){
0563           std::cout << "nz #" << nz << " adc " << adcval[nphi][nz] << " n zero " << zero_count << " last (" << rhitset->m_tpchits[nphi].size() << ") " << ((int)rhitset->m_tpchits[nphi].back()) << std::endl;
0564         if(rhitset->m_tpchits[nphi].size()>=3)
0565         std::cout << " [ " << rhitset->m_tpchits[nphi].size()-1 << ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-1]) << " ] "
0566                     << " [ " << rhitset->m_tpchits[nphi].size()-2 << ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-2]) << " ] "
0567                     << " [ " << rhitset->m_tpchits[nphi].size()-3<< ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-3]) << " ] "
0568                     << std::endl;
0569                     }
0570         */
0571       }
0572       if (zero_count > 0 && (*(rhitset->getHits(nphi))).back() != 0)
0573       {
0574         (*(rhitset->getHits(nphi))).push_back(0);  // pad zero at the end
0575       }
0576       //      std::cout << "sector: " << sector << " lay: " << layer << " nphi: " << nphi << "|" << NPhiBins<< " nzfilled: " << rhitset->m_tpchits[nphi].size() << std::endl;
0577     }
0578 //    count++;
0579 
0580     //    std::cout << "unpacking..." << std::endl;
0581     //    unpack and check for correctness
0582     std::vector<std::vector<uint8_t>> outval;
0583     outval.resize(NPhiBins);
0584     for (int nphi = 0; nphi < NPhiBins; nphi++)
0585     {
0586       outval[nphi].resize(NZBins, 0);
0587       for (int nz = 0; nz < NZBins; nz++)
0588       {
0589         if (outval[nphi][nz] != 0)
0590         {
0591           std::cout << "WARNING!" << std::endl;
0592         }
0593       }
0594     }
0595     //   now we have a clean output array
0596     for (int nphi = 0; nphi < NPhiBins; nphi++)
0597     {
0598       if (rhitset->size(nphi) == 0)
0599       {
0600         continue;
0601       }
0602 
0603       int pindex = 0;
0604       for (unsigned int nzo = 0; nzo < rhitset->size(nphi); nzo++)
0605       {
0606         uint8_t val = (*(rhitset->getHits(nphi)))[nzo];
0607 
0608         if (val == 0)
0609         {
0610           pindex++;
0611         }
0612         else
0613         {
0614           if (nzo == 0)
0615           {
0616             outval[nphi][pindex++] = val;
0617           }
0618           else
0619           {
0620             if (((*(rhitset->getHits(nphi)))[nzo - 1] == 0) && ((*(rhitset->getHits(nphi)))[nzo + 1] == 0))
0621             {  // found zero count
0622               pindex += val;
0623             }
0624             else
0625             {
0626               outval[nphi][pindex++] = val;
0627             }
0628           }
0629         }
0630       }
0631     }
0632     {
0633       // std::cout << "checking..." << std::endl;
0634       for (int nphi = 0; nphi < NPhiBins; nphi++)
0635       {
0636         // std::cout << "nphi: " << nphi << "|" << NPhiBins<< std::endl;
0637         // std::cout << " adcval size "  << adcval[nphi].size() << std::endl;
0638         // std::cout << " outval size "  << outval[nphi].size() << std::endl;
0639         int bad = 0;
0640         for (int nz = 0; nz < NZBins; nz++)
0641         {
0642           ncheck++;
0643           if (adcval[nphi][nz] != outval[nphi][nz])
0644           {
0645             if (outval[nphi][nz] != 255)
0646             {
0647               bad++;
0648               allbad++;
0649             }
0650           }
0651         }
0652 
0653         if (bad > 0)
0654         {
0655           std::cout << "sector: " << sector << " lay: " << layer << " # " << nphi << " bad:" << bad << " packsize: " << (int) rhitset->size(nphi) << std::endl;
0656 
0657           for (int nz = 0; nz < NZBins; nz++)
0658           {
0659             std::cout << "sector: " << sector << " lay: " << layer << " # " << nphi << "|" << nz << " bad:" << bad << " org: " << (int) adcval[nphi][nz] << " pack: " << (int) outval[nphi][nz] << std::endl;
0660           }
0661           for (int nz = 0; nz < (int) rhitset->size(nphi); nz++)
0662           {
0663             std::cout << "lay: " << layer << " # " << nphi << "|" << nz << " bad:" << bad << " pack: " << (int) (*(rhitset->getHits(nphi)))[nz] << std::endl;
0664           }
0665         }
0666       }
0667     }
0668 
0669     //#ifdef WORK     //    if(layer==7)
0670   }
0671   std::cout << " number of hits checked " << ncheck << " bad: " << allbad << std::endl;
0672 
0673   return Fun4AllReturnCodes::EVENT_OK;
0674 }
0675 
0676 int TpcRawWriter::End(PHCompositeNode * /*topNode*/)
0677 {
0678   return Fun4AllReturnCodes::EVENT_OK;
0679 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!