Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:18:03 

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/PHG4TpcCylinderGeom.h>
0026 #include <g4detectors/PHG4TpcCylinderGeomContainer.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   PHG4TpcCylinderGeomContainer *geom_container =
0185       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0186   if (!geom_container)
0187   {
0188     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << 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   // loop over the TPC HitSet objects
0315   TrkrHitSetContainer::ConstRange tpc_hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0316   //  const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
0317 
0318   int ncheck = 0;
0319   int allbad = 0;
0320 
0321   for (TrkrHitSetContainer::ConstIterator hitsetitr = tpc_hitsetrange.first;
0322        hitsetitr != tpc_hitsetrange.second;
0323        ++hitsetitr)
0324   // TrkrHitSetContainer::ConstIterator hitsetitr = tpc_hitsetrange.first;
0325   {
0326     TrkrHitSet *hitset = hitsetitr->second;
0327     unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0328     //   if(layer!=7) continue;
0329     //    int side = TpcDefs::getSide(hitsetitr->first);
0330     unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0331     PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
0332     unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
0333     unsigned short NPhiBinsSector = NPhiBins / 12;
0334     unsigned short NZBins = (unsigned short) layergeom->get_zbins();
0335     // unsigned short NZBinsSide = NZBins/2;
0336     unsigned short NZBinsSide = NZBins;
0337     unsigned short NZBinsMin = 0;
0338     unsigned short PhiOffset = NPhiBinsSector * sector;
0339     unsigned short ZOffset = NZBinsMin;
0340     std::vector<int> nhits_thispad;
0341     std::vector<std::vector<unsigned short>> adcval;
0342     adcval.resize(NPhiBins);
0343     for (int nphi = 0; nphi < NPhiBins; nphi++)
0344     {
0345       for (int nz = 0; nz < NZBins; nz++)
0346       {
0347         adcval[nphi].push_back(0);
0348         nhits_thispad.push_back(0);
0349       }
0350     }
0351     // std::vector<std::vector<unsigned short>> outadc;
0352     // outadc.resize(NPhiBins);
0353 
0354     //    RawHitSetContainerv1::Iterator rawhitset_iter =
0355     m_rawhits->findOrAddHitSet(hitsetitr->first);
0356     RawHitSetv1 *rhitset = dynamic_cast<RawHitSetv1 *>(m_rawhits->findHitSet(hitsetitr->first));
0357     rhitset->setTpcPhiBins(NPhiBins);
0358 
0359     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0360     int zbinmax = 498;
0361     int zbinmin = 0;
0362     if (layer >= 7 && layer < 22)
0363     {
0364       int etacut = 249 - ((50 + (layer - 7)) / 105.5) * 249;
0365       zbinmin = etacut;
0366       zbinmax -= etacut;
0367     }
0368     if (layer >= 22 && layer <= 48)
0369     {
0370       int etacut = 249 - ((65 + ((40.5 / 26) * (layer - 22))) / 105.5) * 249;
0371       zbinmin = etacut;
0372       zbinmax -= etacut;
0373     }
0374 
0375     // std::cout << " layer: " << layer << " zbin limit " << zbinmin << " | " << zbinmax <<std::endl;
0376     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0377          hitr != hitrangei.second;
0378          ++hitr)
0379     {
0380       if (TpcDefs::getPad(hitr->first) - PhiOffset < 0)
0381       {
0382         continue;
0383       }
0384       if (TpcDefs::getTBin(hitr->first) - ZOffset < 0)
0385       {
0386         // std::cout << "WARNING zbin out of range: " << TpcDefs::getTBin(hitr->first) - zoffset  << " | " << zbins <<std::endl;
0387         continue;
0388       }
0389       unsigned short phibin = TpcDefs::getPad(hitr->first) - PhiOffset;
0390       unsigned short zbin = TpcDefs::getTBin(hitr->first) - ZOffset;
0391       unsigned short zbinorg = TpcDefs::getTBin(hitr->first);
0392       if (phibin >= NPhiBinsSector)
0393       {
0394         // std::cout << "WARNING phibin out of range: " << phibin << " | " << phibins << std::endl;
0395         continue;
0396       }
0397       if (zbin >= NZBinsSide)
0398       {
0399         // std::cout << "WARNING z bin out of range: " << zbin << " | " << zbins << std::endl;
0400         continue;
0401       }
0402       if (zbinorg > zbinmax || zbinorg < zbinmin)
0403       {
0404         continue;
0405       }
0406       float_t fadc = (hitr->second->getAdc()) - pedestal;  // proper int rounding +0.5
0407       unsigned short adc = 0;
0408       if (fadc > 0)
0409       {
0410         adc = (unsigned short) fadc;
0411       }
0412 
0413       if (adc > 0)
0414       {
0415         nhits_thispad[phibin]++;
0416         unsigned short it = TpcDefs::getTBin(hitr->first);
0417         adcval[phibin][it] = adc;
0418         // if(layer==20&&sector==4&&phibin==77)
0419         // std::cout << "input sector: " << sector << " lay: " << layer << " # " << phibin << "|"<< it <<" adc: " << (int)adcval[phibin][it] << " pack: " << std::endl;
0420         /*RawHitTpc  *rHit = new RawHitTpc;
0421         unsigned short it = TpcDefs::getTBin(hitr->first);
0422         rHit->setTBin(it);
0423         unsigned short iadc = adc;
0424         rHit->setAdc(iadc);
0425         rhitset->addTpcHit(phibin,rHit);
0426         */
0427       }
0428     }
0429     //    assert(adcval);
0430     for (int nphi = 0; nphi < NPhiBins; nphi++)
0431     {
0432       for (int nz = 0; nz < NZBins; nz++)
0433       {
0434         /*if(nphi>=1&&nz<NZBins-1)//Remove isolated single timebin hits
0435           if(layer==29&&sector==0&&nphi==0){
0436             if(nz>209&&nz<215){
0437               std::cout << "single check sector: " << sector << " lay: " << layer << " # " << nphi << "|"<< nz <<" adc: " << (int)adcval[nphi][nz] << " pack: " << std::endl;
0438             }
0439           }
0440         */
0441         if (adcval[nphi][nz] > 0)
0442         {
0443           if ((adcval[nphi][nz - 1] == 0) && (adcval[nphi][nz + 1] == 0))
0444           {
0445             // std::cout << "single hit sector: " << sector << " lay: " << layer << " # " << nphi << "|"<< nz <<" adc: " << (int)adcval[nphi][nz] << " pack: " << std::endl;
0446             adcval[nphi][nz] = 0;
0447           }
0448         }
0449       }
0450     }
0451     for (int nphi = 0; nphi < NPhiBins; nphi++)
0452     {
0453       int zero_count = 0;
0454       // int outpos = 0;
0455       if (nhits_thispad[nphi] == 0)
0456       {
0457         continue;
0458       }
0459       for (int nz = 0; nz < NZBins; nz++)
0460       {
0461         uint8_t thisadc = 0;
0462         if (nphi >= 1 && nz < NZBins - 1)
0463         {  // Remove isolated single timebin hits
0464           if (adcval[nphi][nz] > 0)
0465           {
0466             if ((adcval[nphi][nz - 1] == 0) && (adcval[nphi][nz + 1] == 0))
0467             {
0468               adcval[nphi][nz] = 0;
0469             }
0470           }
0471         }
0472 
0473         if (adcval[nphi][nz] > 255)
0474         {  // take care of overflows, we want to store 8 bit ADC values
0475           thisadc = 255;
0476           //        adcval[nphi][nz]=255;
0477         }
0478         else
0479         {
0480           thisadc = adcval[nphi][nz];
0481         }
0482         // std::cout << "nz: " << nz << " adc " << (int)thisadc << std::endl;
0483         if (thisadc > 0)
0484         {  // non zero adc, fill adc
0485           zero_count = 0;
0486           (*(rhitset->getHits(nphi))).push_back(thisadc);
0487           if (layer == 222 && sector == 6 && nphi == 10)
0488           {
0489             std::cout << "0#" << nphi << "nz= " << nz << " filling " << (int) thisadc << " at " << rhitset->size(nphi) << std::endl;
0490           }
0491         }
0492         else
0493         {
0494           zero_count++;
0495           //      if(nz==0){
0496           //  rhitset->m_tpchits[nphi].push_back(0);
0497           // std::cout << " 1filling " << 0 << "|" << rhitset->m_tpchits[nphi].back() << " | " << rhitset->m_tpchits[nphi][outpos++] << std::endl;
0498           // if(layer==222&&sector==6&&nphi==10)std::cout << "1#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
0499           // }else
0500           /*
0501           if(adcval[nphi][nz-1]>0){ //trailing edge 0
0502             rhitset->m_tpchits[nphi].push_back(0);
0503             if(layer==222&&sector==6&&nphi==10)std::cout << "2#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
0504             //std::cout << " 2filling " << 0 << "|" << rhitset->m_tpchits[nphi].back() << " | " << rhitset->m_tpchits[nphi][outpos++] << std::endl;
0505           }
0506           */
0507         }
0508         if (zero_count > 1)
0509         {
0510           if (nz < NZBins - 2)
0511           {
0512             if (adcval[nphi][nz + 2] > 0 && adcval[nphi][nz + 1] == 0)
0513             {  // fill zero count, end of zero series
0514               (*(rhitset->getHits(nphi))).push_back(zero_count - 1);
0515               zero_count = 0;
0516               if (layer == 222 && sector == 6 && nphi == 10)
0517               {
0518                 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;
0519               }
0520             }
0521             if (adcval[nphi][nz + 1] > 0)
0522             {
0523               (*(rhitset->getHits(nphi))).push_back(0);  // leading edge zero
0524               if (layer == 222 && sector == 6 && nphi == 10)
0525               {
0526                 std::cout << "4#" << nphi << "nz= " << nz << " filling " << 0 << " at " << rhitset->size(nphi) << std::endl;
0527               }
0528 
0529               // std::cout << " 4filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
0530             }
0531           }
0532         }
0533         else
0534         {
0535           if (zero_count == 1)
0536           {
0537             (*(rhitset->getHits(nphi))).push_back(0);
0538             if (layer == 222 && sector == 6 && nphi == 10)
0539             {
0540               std::cout << "5#" << nphi << "nz= " << nz << " filling " << 0 << " at " << rhitset->size(nphi) << std::endl;
0541             }
0542             // std::cout << " 5filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
0543           }
0544         }
0545         if (zero_count == 254)
0546         {
0547           zero_count = 0;
0548           (*(rhitset->getHits(nphi))).push_back(254 - 2);
0549           (*(rhitset->getHits(nphi))).push_back(0);
0550           if (layer == 222 && sector == 6 && nphi == 10)
0551           {
0552             std::cout << "6#" << nphi << "nz= " << nz << " filling " << 254 - 1 << " at " << rhitset->size(nphi) << std::endl;
0553           }
0554           // std::cout << "6 filling " << 254-1 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
0555           //  rhitset->m_tpchits[nphi].push_back(0);
0556           //  if(layer==222&&sector==6&&nphi==10)std::cout << "7#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
0557           // std::cout << " 7filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
0558         }
0559         /*
0560           if(layer%15==0){
0561           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;
0562         if(rhitset->m_tpchits[nphi].size()>=3)
0563         std::cout << " [ " << rhitset->m_tpchits[nphi].size()-1 << ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-1]) << " ] "
0564                     << " [ " << rhitset->m_tpchits[nphi].size()-2 << ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-2]) << " ] "
0565                     << " [ " << rhitset->m_tpchits[nphi].size()-3<< ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-3]) << " ] "
0566                     << std::endl;
0567                     }
0568         */
0569       }
0570       if (zero_count > 0 && (*(rhitset->getHits(nphi))).back() != 0)
0571       {
0572         (*(rhitset->getHits(nphi))).push_back(0);  // pad zero at the end
0573       }
0574       //      std::cout << "sector: " << sector << " lay: " << layer << " nphi: " << nphi << "|" << NPhiBins<< " nzfilled: " << rhitset->m_tpchits[nphi].size() << std::endl;
0575     }
0576 //    count++;
0577 
0578     //    std::cout << "unpacking..." << std::endl;
0579     //    unpack and check for correctness
0580     std::vector<std::vector<uint8_t>> outval;
0581     outval.resize(NPhiBins);
0582     for (int nphi = 0; nphi < NPhiBins; nphi++)
0583     {
0584       outval[nphi].resize(NZBins, 0);
0585       for (int nz = 0; nz < NZBins; nz++)
0586       {
0587         if (outval[nphi][nz] != 0)
0588         {
0589           std::cout << "WARNING!" << std::endl;
0590         }
0591       }
0592     }
0593     //   now we have a clean output array
0594     for (int nphi = 0; nphi < NPhiBins; nphi++)
0595     {
0596       if (rhitset->size(nphi) == 0)
0597       {
0598         continue;
0599       }
0600 
0601       int pindex = 0;
0602       for (unsigned int nzo = 0; nzo < rhitset->size(nphi); nzo++)
0603       {
0604         uint8_t val = (*(rhitset->getHits(nphi)))[nzo];
0605 
0606         if (val == 0)
0607         {
0608           pindex++;
0609         }
0610         else
0611         {
0612           if (nzo == 0)
0613           {
0614             outval[nphi][pindex++] = val;
0615           }
0616           else
0617           {
0618             if (((*(rhitset->getHits(nphi)))[nzo - 1] == 0) && ((*(rhitset->getHits(nphi)))[nzo + 1] == 0))
0619             {  // found zero count
0620               pindex += val;
0621             }
0622             else
0623             {
0624               outval[nphi][pindex++] = val;
0625             }
0626           }
0627         }
0628       }
0629     }
0630     {
0631       // std::cout << "checking..." << std::endl;
0632       for (int nphi = 0; nphi < NPhiBins; nphi++)
0633       {
0634         // std::cout << "nphi: " << nphi << "|" << NPhiBins<< std::endl;
0635         // std::cout << " adcval size "  << adcval[nphi].size() << std::endl;
0636         // std::cout << " outval size "  << outval[nphi].size() << std::endl;
0637         int bad = 0;
0638         for (int nz = 0; nz < NZBins; nz++)
0639         {
0640           ncheck++;
0641           if (adcval[nphi][nz] != outval[nphi][nz])
0642           {
0643             if (outval[nphi][nz] != 255)
0644             {
0645               bad++;
0646               allbad++;
0647             }
0648           }
0649         }
0650 
0651         if (bad > 0)
0652         {
0653           std::cout << "sector: " << sector << " lay: " << layer << " # " << nphi << " bad:" << bad << " packsize: " << (int) rhitset->size(nphi) << std::endl;
0654 
0655           for (int nz = 0; nz < NZBins; nz++)
0656           {
0657             std::cout << "sector: " << sector << " lay: " << layer << " # " << nphi << "|" << nz << " bad:" << bad << " org: " << (int) adcval[nphi][nz] << " pack: " << (int) outval[nphi][nz] << std::endl;
0658           }
0659           for (int nz = 0; nz < (int) rhitset->size(nphi); nz++)
0660           {
0661             std::cout << "lay: " << layer << " # " << nphi << "|" << nz << " bad:" << bad << " pack: " << (int) (*(rhitset->getHits(nphi)))[nz] << std::endl;
0662           }
0663         }
0664       }
0665     }
0666 
0667     //#ifdef WORK     //    if(layer==7)
0668   }
0669   std::cout << " number of hits checked " << ncheck << " bad: " << allbad << std::endl;
0670 
0671   return Fun4AllReturnCodes::EVENT_OK;
0672 }
0673 
0674 int TpcRawWriter::End(PHCompositeNode * /*topNode*/)
0675 {
0676   return Fun4AllReturnCodes::EVENT_OK;
0677 }
[ 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!