Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "TpcCombinedRawDataUnpacker.h"
0002 
0003 #include <trackbase/TpcDefs.h>
0004 #include <trackbase/TrkrDefs.h>  // for hitkey, hitsetkey
0005 #include <trackbase/TrkrHit.h>
0006 #include <trackbase/TrkrHitSet.h>
0007 #include <trackbase/TrkrHitSetContainer.h>
0008 #include <trackbase/TrkrHitSetContainerv1.h>
0009 #include <trackbase/TrkrHitv2.h>
0010 
0011 #include <g4detectors/PHG4TpcCylinderGeom.h>
0012 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0013 
0014 #include <ffarawobjects/TpcRawHit.h>
0015 #include <ffarawobjects/TpcRawHitContainer.h>
0016 
0017 #include <cdbobjects/CDBTTree.h>
0018 
0019 #include <ffamodules/CDBInterface.h>
0020 
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/Fun4AllServer.h>
0023 
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/PHIODataNode.h>  // for PHIODataNode
0026 #include <phool/PHNodeIterator.h>
0027 #include <phool/PHObject.h>  // for PHObject
0028 #include <phool/getClass.h>
0029 #include <phool/phool.h>  // for PHWHERE
0030 
0031 #include <TFile.h>
0032 #include <TH1.h>
0033 #include <TH2.h>
0034 #include <TNtuple.h>
0035 #include <TSystem.h>
0036 
0037 #include <algorithm>
0038 #include <cmath>
0039 #include <cstdint>   // for exit
0040 #include <cstdlib>   // for exit
0041 #include <iostream>  // for operator<<, endl, bas...
0042 #include <map>       // for _Rb_tree_iterator
0043 #include <memory>
0044 #include <utility>
0045 
0046 TpcCombinedRawDataUnpacker::TpcCombinedRawDataUnpacker(std::string const& name, std::string const& outF)
0047   : SubsysReco(name)
0048   , outfile_name(outF)
0049 {
0050   // Do nothing
0051 }
0052 TpcCombinedRawDataUnpacker::~TpcCombinedRawDataUnpacker()
0053 {
0054   delete m_cdbttree;
0055 }
0056 void TpcCombinedRawDataUnpacker::ReadZeroSuppressedData()
0057 {
0058   m_do_zs_emulation = true;
0059   m_do_baseline_corr = false;
0060   auto cdb = CDBInterface::instance();
0061   std::string dir = cdb->getUrl("TPC_ZS_THRESHOLDS");
0062 
0063   auto cdbtree = std::make_unique<CDBTTree>(dir);
0064   std::ostringstream name;
0065   if(dir.empty())
0066  {
0067   if(Verbosity() > 1)
0068   {
0069     std::cout << "use default tpc zs threshold of 20" << std::endl;
0070   }
0071   return;
0072  }
0073  
0074   for(int i=0; i<3; i++)
0075   {
0076     name.str("");
0077     name << "R"<<i+1<<"ADUThreshold";
0078     m_zs_threshold[i] = cdbtree->GetSingleFloatValue(name.str().c_str());
0079     if(Verbosity() > 1)
0080     {
0081       std::cout << "Loading ADU threshold of " << m_zs_threshold[i] << " for region " << i << std::endl;
0082     }
0083   }
0084 
0085 }
0086 int TpcCombinedRawDataUnpacker::Init(PHCompositeNode* /*topNode*/)
0087 {
0088   std::cout << "TpcCombinedRawDataUnpacker::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0089 
0090   m_cdb = CDBInterface::instance();
0091   std::string calibdir = m_cdb->getUrl("TPC_FEE_CHANNEL_MAP");
0092 
0093   if (calibdir[0] == '/')
0094   {
0095     // use generic CDBTree to load
0096     m_cdbttree = new CDBTTree(calibdir);
0097     m_cdbttree->LoadCalibrations();
0098   }
0099   else
0100   {
0101     std::cout << "TpcRawDataDecoder::::InitRun No calibration file found" << std::endl;
0102     exit(1);
0103   }
0104 
0105   return Fun4AllReturnCodes::EVENT_OK;
0106 }
0107 
0108 int TpcCombinedRawDataUnpacker::InitRun(PHCompositeNode* topNode)
0109 {
0110   if (!topNode)
0111   {
0112     std::cout << "TpcCombinedRawDataUnpacker::InitRun(PHCompositeNode* topNode)" << std::endl;
0113     std::cout << "\tCould not retrieve topNode; doing nothing" << std::endl;
0114     exit(1);
0115     gSystem->Exit(1);
0116 
0117     return 1;
0118   }
0119 
0120   PHNodeIterator dst_itr(topNode);
0121   PHCompositeNode* dst_node = dynamic_cast<PHCompositeNode*>(dst_itr.findFirst("PHCompositeNode", "DST"));
0122   if (!dst_node)
0123   {
0124     if (Verbosity())
0125     {
0126       std::cout << "TpcCombinedRawDataUnpacker::InitRun(PHCompositeNode* topNode)" << std::endl;
0127     }
0128     if (Verbosity())
0129     {
0130       std::cout << "\tCould not retrieve dst_node; doing nothing" << std::endl;
0131     }
0132     exit(1);
0133     gSystem->Exit(1);
0134 
0135     return 1;
0136   }
0137 
0138   PHNodeIterator trkr_itr(dst_node);
0139   PHCompositeNode* trkr_node = dynamic_cast<PHCompositeNode*>(trkr_itr.findFirst("PHCompositeNode", "TRKR"));
0140   if (!trkr_node)
0141   {
0142     trkr_node = new PHCompositeNode("TRKR");
0143     dst_node->addNode(trkr_node);
0144   }
0145 
0146   TrkrHitSetContainer* trkr_hit_set_container = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0147   if (!trkr_hit_set_container)
0148   {
0149     if (Verbosity())
0150     {
0151       std::cout << "TpcCombinedRawDataUnpacker::InitRun(PHCompositeNode* topNode)" << std::endl;
0152     }
0153     if (Verbosity())
0154     {
0155       std::cout << "\tMaking TrkrHitSetContainer" << std::endl;
0156     }
0157 
0158     trkr_hit_set_container = new TrkrHitSetContainerv1;
0159     PHIODataNode<PHObject>* new_node = new PHIODataNode<PHObject>(trkr_hit_set_container, "TRKR_HITSET", "PHObject");
0160     trkr_node->addNode(new_node);
0161   }
0162 
0163   TpcRawHitContainer* tpccont = findNode::getClass<TpcRawHitContainer>(topNode, m_TpcRawNodeName);
0164   if (!tpccont)
0165   {
0166     std::cout << PHWHERE << std::endl;
0167     std::cout << "TpcCombinedRawDataUnpacker::process_event(PHCompositeNode* topNode)" << std::endl;
0168     std::cout << "Could not get \"" << m_TpcRawNodeName << "\" from Node Tree" << std::endl;
0169     std::cout << "Removing module" << std::endl;
0170 
0171     Fun4AllServer* se = Fun4AllServer::instance();
0172     se->unregisterSubsystem(this);
0173     return Fun4AllReturnCodes::EVENT_OK;
0174   }
0175 
0176   if (m_writeTree)
0177   {
0178     m_file = new TFile(outfile_name.c_str(), "RECREATE");
0179     m_ntup = new TNtuple("NT", "NT", "event:gtmbco:packid:ep:sector:side:fee:rx:entries:ped:width");
0180     m_ntup_hits = new TNtuple("NTH", "NTH", "event:gtmbco:packid:ep:sector:side:fee:chan:sampadd:sampch:phibin:tbin:layer:adc:ped:width");
0181     m_ntup_hits_corr = new TNtuple("NTC", "NTC", "event:gtmbco:packid:ep:sector:side:fee:chan:sampadd:sampch:phibin:tbin:layer:adc:ped:width:corr");
0182     if (m_doChanHitsCut)
0183     {
0184       m_HitChanDis = new TH2F("HitChanDis", "HitChanDis", 451, -0.5, 450.5, 256, -0.5, 255.5);
0185       m_HitsinChan = new TH1F("HitsinChan", "HitsinChan", 451, -0.5, 450.5);
0186     }
0187   }
0188 
0189   if (Verbosity() >= 1)
0190   {
0191     std::cout << "TpcCombinedRawDataUnpacker:: _ped_sig_cut = " << m_ped_sig_cut << std::endl;
0192     std::cout << "TpcCombinedRawDataUnpacker:: startevt = " << startevt << std::endl;
0193     std::cout << "TpcCombinedRawDataUnpacker:: endevt = " << endevt << std::endl;
0194   }
0195 
0196   // check run number if presamples need to be shifted, which went from 80 -> 120
0197   // at 41624
0198   Fun4AllServer* se = Fun4AllServer::instance();
0199   if (se->RunNumber() < 41624)
0200   {
0201     m_presampleShift = 0;
0202   }
0203 
0204   return Fun4AllReturnCodes::EVENT_OK;
0205 }
0206 
0207 int TpcCombinedRawDataUnpacker::process_event(PHCompositeNode* topNode)
0208 {
0209   if (_ievent < startevt || _ievent > endevt)
0210   {
0211     if (Verbosity() > 1)
0212     {
0213       std::cout << " Skip event " << _ievent << std::endl;
0214     }
0215     _ievent++;
0216     return Fun4AllReturnCodes::DISCARDEVENT;
0217   }
0218   _ievent++;
0219 
0220   TrkrHitSetContainer* trkr_hit_set_container = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0221   if (!trkr_hit_set_container)
0222   {
0223     std::cout << PHWHERE << std::endl;
0224     std::cout << "TpcCombinedRawDataUnpacker::process_event(PHCompositeNode* topNode)" << std::endl;
0225     std::cout << "Could not get \"TRKR_HITSET\" from Node Tree" << std::endl;
0226     std::cout << "Exiting" << std::endl;
0227     gSystem->Exit(1);
0228     exit(1);
0229 
0230     return Fun4AllReturnCodes::DISCARDEVENT;
0231   }
0232 
0233   TpcRawHitContainer* tpccont = findNode::getClass<TpcRawHitContainer>(topNode, m_TpcRawNodeName);
0234   if (!tpccont)
0235   {
0236     std::cout << PHWHERE << std::endl;
0237     std::cout << "TpcCombinedRawDataUnpacker::process_event(PHCompositeNode* topNode)" << std::endl;
0238     std::cout << "Could not get \"" << m_TpcRawNodeName << "\" from Node Tree" << std::endl;
0239     std::cout << "Exiting" << std::endl;
0240 
0241     gSystem->Exit(1);
0242     exit(1);
0243   }
0244 
0245   PHG4TpcCylinderGeomContainer* geom_container =
0246       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0247   if (!geom_container)
0248   {
0249     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0250     return Fun4AllReturnCodes::ABORTRUN;
0251   }
0252 
0253   TrkrDefs::hitsetkey hit_set_key = 0;
0254   TrkrDefs::hitkey hit_key = 0;
0255   TrkrHitSetContainer::Iterator hit_set_container_itr;
0256   TrkrHit* hit = nullptr;
0257 
0258   uint64_t bco_min = UINT64_MAX;
0259   uint64_t bco_max = 0;
0260 
0261   const auto nhits = tpccont->get_nhits();
0262 
0263   int max_time_range = 0;
0264 
0265   for (unsigned int i = 0; i < nhits; i++)
0266   {
0267     TpcRawHit* tpchit = tpccont->get_hit(i);
0268     uint64_t gtm_bco = tpchit->get_gtm_bco();
0269 
0270     bco_min = std::min(gtm_bco, bco_min);
0271     bco_max = std::max(gtm_bco, bco_max);
0272 
0273     int fee = tpchit->get_fee();
0274     int channel = tpchit->get_channel();
0275 
0276     int feeM = FEE_map[fee];
0277     if (FEE_R[fee] == 2)
0278     {
0279       feeM += 6;
0280     }
0281     if (FEE_R[fee] == 3)
0282     {
0283       feeM += 14;
0284     }
0285 
0286     int side = 1;
0287     int32_t packet_id = tpchit->get_packetid();
0288     int ep = (packet_id - 4000) % 10;
0289     int sector = (packet_id - 4000 - ep) / 10;
0290     if (sector > 11)
0291     {
0292       side = 0;
0293     }
0294 
0295     unsigned int key = (256 * (feeM)) + channel;
0296     std::string varname = "layer";
0297     int layer = m_cdbttree->GetIntValue(key, varname);
0298     // antenna pads will be in 0 layer
0299     if (layer <= 6)
0300     {
0301       continue;
0302     }
0303 
0304     uint16_t sampadd = tpchit->get_sampaaddress();
0305     uint16_t sampch = tpchit->get_sampachannel();
0306     //    uint16_t sam = tpchit->get_samples();
0307     max_time_range = tpchit->get_samples();
0308     varname = "phi";  // + std::to_string(key);
0309     double phi = ((side == 1 ? 1 : -1) * (m_cdbttree->GetDoubleValue(key, varname) - M_PI / 2.)) + ((sector % 12) * M_PI / 6);
0310     PHG4TpcCylinderGeom* layergeom = geom_container->GetLayerCellGeom(layer);
0311     unsigned int phibin = layergeom->get_phibin(phi, side);
0312     unsigned int region = 0;
0313     if(layer > 15)
0314     {
0315       region = 1;
0316     }
0317     if( layer > 31)
0318     {
0319       region = 2;
0320     }
0321     hit_set_key = TpcDefs::genHitSetKey(layer, (mc_sectors[sector % 12]), side);
0322     hit_set_container_itr = trkr_hit_set_container->findOrAddHitSet(hit_set_key);
0323 
0324     float hpedestal = 0;
0325     float hpedwidth = 0;
0326 
0327     if (Verbosity() > 2)
0328     {
0329       std::cout << "TpcCombinedRawDataUnpacker:: do zero suppression" << std::endl;
0330     }
0331     TH2* feehist = nullptr;
0332     hpedestal = 60;
0333     hpedwidth = m_zs_threshold[region];
0334 
0335     unsigned int pad_key = create_pad_key(side, layer, phibin);
0336 
0337     std::map<unsigned int, chan_info>::iterator chan_it = chan_map.find(pad_key);
0338     if (chan_it != chan_map.end())
0339     {
0340       (*chan_it).second.ped = hpedestal;
0341       (*chan_it).second.width = hpedwidth;
0342     }
0343     else
0344     {
0345       chan_info nucinfo;
0346       nucinfo.fee = fee;
0347       nucinfo.ped = hpedestal;
0348       nucinfo.width = hpedwidth;
0349       chan_map.insert(std::make_pair(pad_key, nucinfo));
0350     }
0351     int rx = get_rx(layer);
0352     unsigned int fee_key = create_fee_key(side, mc_sectors[sector % 12], rx, fee);
0353     // find or insert TH2C;
0354     std::map<unsigned int, TH2*>::iterator fee_map_it;
0355 
0356     fee_map_it = feeadc_map.find(fee_key);
0357     if (fee_map_it != feeadc_map.end())
0358     {
0359       feehist = (*fee_map_it).second;
0360     }
0361     else
0362     {
0363       std::string histname = "h" + std::to_string(fee_key);
0364       feehist = new TH2C(histname.c_str(), "histname", max_time_range + 1, -0.5, max_time_range + 0.5, 501, -0.5, 1000.5);
0365       feeadc_map.insert(std::make_pair(fee_key, feehist));
0366       std::vector<int> feeentries(feehist->GetNbinsX(), 0);
0367       feeentries_map.insert(std::make_pair(fee_key, feeentries));
0368     }
0369     auto fee_entries_it = feeentries_map.find(fee_key);
0370     std::vector<int>& fee_entries_vec = (*fee_entries_it).second;
0371 
0372     float threshold_cut = m_zs_threshold[region];
0373 
0374     int nhitschan = 0;
0375 
0376     if (m_doChanHitsCut)
0377     {
0378       for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(tpchit->CreateAdcIterator()); !adc_iterator->IsDone(); adc_iterator->Next())
0379       {
0380         const uint16_t s = adc_iterator->CurrentTimeBin();
0381         const uint16_t adc = adc_iterator->CurrentAdc();
0382         int t = s - m_presampleShift - m_t0;
0383         if (t < 0)
0384         {
0385           continue;
0386         }
0387         if (feehist != nullptr)
0388         {
0389           if (adc > 0)
0390           {
0391             if ((float(adc) - hpedestal) > threshold_cut)
0392             {
0393               nhitschan++;
0394             }
0395           }
0396         }
0397       }
0398       if (m_writeTree)
0399       {
0400         m_HitChanDis->Fill(nhitschan, channel);
0401         m_HitsinChan->Fill(nhitschan);
0402       }
0403       if (nhitschan > m_ChanHitsCut)
0404       {
0405         continue;
0406       }
0407     }
0408 
0409     for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(tpchit->CreateAdcIterator());
0410          !adc_iterator->IsDone();
0411          adc_iterator->Next())
0412     {
0413       const uint16_t s = adc_iterator->CurrentTimeBin();
0414       const uint16_t adc = adc_iterator->CurrentAdc();
0415       int t = s - m_presampleShift - m_t0;
0416       if (t < 0)
0417       {
0418         continue;
0419       }
0420       if (feehist != nullptr)
0421       {
0422         if (adc > 0)
0423         {
0424           if ((float(adc) - hpedestal) > threshold_cut)
0425           {
0426             feehist->Fill(t, adc - hpedestal);
0427             if (t < (int) fee_entries_vec.size())
0428             {
0429               fee_entries_vec[t]++;
0430             }
0431           }
0432         }
0433       }
0434 
0435       if ((float(adc) - hpedestal) > threshold_cut)
0436       {
0437         hit_key = TpcDefs::genHitKey(phibin, (unsigned int) t);
0438         // find existing hit, or create new one
0439         hit = hit_set_container_itr->second->getHit(hit_key);
0440         if (!hit)
0441         {
0442           hit = new TrkrHitv2();
0443           hit->setAdc(float(adc) - hpedestal);
0444           hit_set_container_itr->second->addHitSpecificKey(hit_key, hit);
0445         }
0446 
0447         if (m_writeTree)
0448         {
0449           float fXh[18];
0450           int nh = 0;
0451 
0452           fXh[nh++] = _ievent - 1;
0453           fXh[nh++] = gtm_bco;                  // gtm_bco;
0454           fXh[nh++] = packet_id;                // packet_id;
0455           fXh[nh++] = ep;                       // ep;
0456           fXh[nh++] = mc_sectors[sector % 12];  // Sector;
0457           fXh[nh++] = side;
0458           fXh[nh++] = fee;
0459           fXh[nh++] = channel;  // channel;
0460           fXh[nh++] = sampadd;  // sampadd;
0461           fXh[nh++] = sampch;   // sampch;
0462           fXh[nh++] = (float) phibin;
0463           fXh[nh++] = (float) t;
0464           fXh[nh++] = layer;
0465           fXh[nh++] = (float(adc) - hpedestal);
0466           fXh[nh++] = hpedestal;
0467           fXh[nh++] = hpedwidth;
0468           m_ntup_hits->Fill(fXh);
0469         }
0470       }
0471     }
0472   }
0473 
0474   if (m_do_baseline_corr == true)
0475   {
0476     // Histos filled now process them for fee local baselines
0477 
0478     int nhistfilled = 0;
0479     int nhisttotal = 0;
0480     for (auto& hiter : feeadc_map)
0481     {
0482       if (hiter.second != nullptr)
0483       {
0484         unsigned int fee_key = hiter.first;
0485         unsigned int side;
0486         unsigned int sector;
0487         unsigned int rx;
0488         unsigned int fee;
0489         unpack_fee_key(side, sector, rx, fee, fee_key);
0490         TH2* hist2d = hiter.second;
0491         std::map<unsigned int, std::vector<int>>::iterator fee_entries_it = feeentries_map.find(fee_key);
0492         if (fee_entries_it == feeentries_map.end())
0493         {
0494           continue;
0495           //      fee_entries_vec_it = (*fee_entries_it).second.begin();
0496         }
0497         std::vector<int>::iterator fee_entries_vec_it = (*fee_entries_it).second.begin();
0498 
0499         std::vector<float> pedvec(hist2d->GetNbinsX(), 0);
0500         feebaseline_map.insert(std::make_pair(hiter.first, pedvec));
0501         std::map<unsigned int, std::vector<float>>::iterator fee_blm_it = feebaseline_map.find(hiter.first);
0502         (*fee_blm_it).second.resize(hist2d->GetNbinsX(), 0);
0503         for (int binx = 1; binx < hist2d->GetNbinsX(); binx++)
0504         {
0505           double timebin = (hist2d->GetXaxis())->GetBinCenter(binx);
0506           std::string histname1d = "h" + std::to_string(hiter.first) + "_" + std::to_string((int) timebin);
0507           nhisttotal++;
0508           float local_ped = 0;
0509           float local_width = 0;
0510           float entries = fee_entries_vec_it[timebin];
0511           if (fee_entries_vec_it[timebin] > 100)
0512           {
0513             nhistfilled++;
0514 
0515             TH1D* hist1d = hist2d->ProjectionY(histname1d.c_str(), binx, binx);
0516             if (hist1d->GetEntries() != fee_entries_vec_it[timebin])
0517             {
0518               std::cout << " vec " << fee_entries_vec_it[timebin]
0519                         << " hist " << hist1d->GetEntries()
0520                         << std::endl;
0521             }
0522             if (hist1d->GetEntries() > 10)
0523             {
0524               int maxbin = hist1d->GetMaximumBin();
0525               // calc peak position
0526               double hadc_sum = 0.0;
0527               double hibin_sum = 0.0;
0528               double hibin2_sum = 0.0;
0529 
0530               for (int isum = -3; isum <= 3; isum++)
0531               {
0532                 float val = hist1d->GetBinContent(maxbin + isum);
0533                 float center = hist1d->GetBinCenter(maxbin + isum);
0534                 hibin_sum += center * val;
0535                 hibin2_sum += center * center * val;
0536                 hadc_sum += val;
0537               }
0538               local_ped = hibin_sum / hadc_sum;
0539               local_width = sqrt((hibin2_sum / hadc_sum) - (local_ped * local_ped));
0540             }
0541             delete hist1d;
0542           }
0543           (*fee_blm_it).second[(int) timebin] = local_ped + m_baseline_nsigma * local_width;
0544 
0545           if (m_writeTree)
0546           {
0547             float fXh[11];
0548             int nh = 0;
0549 
0550             fXh[nh++] = _ievent - 1;
0551             fXh[nh++] = 0;                        // gtm_bco;
0552             fXh[nh++] = 0;                        // packet_id;
0553             fXh[nh++] = 0;                        // ep;
0554             fXh[nh++] = mc_sectors[sector % 12];  // Sector;
0555             fXh[nh++] = side;
0556             fXh[nh++] = fee;
0557             fXh[nh++] = rx;
0558             fXh[nh++] = entries;
0559             fXh[nh++] = local_ped;
0560             fXh[nh++] = local_width;
0561             m_ntup->Fill(fXh);
0562           }
0563         }
0564       }
0565     }
0566     if (Verbosity() >= 1)
0567     {
0568       std::cout << " filled " << nhistfilled
0569                 << " total " << nhisttotal
0570                 << std::endl;
0571 
0572       std::cout << "second loop " << m_do_baseline_corr << std::endl;
0573     }
0574 
0575     // second loop over hits to apply baseline correction
0576     TrkrHitSetContainer::ConstRange hitsetrange;
0577     hitsetrange = trkr_hit_set_container->getHitSets(TrkrDefs::TrkrId::tpcId);
0578 
0579     for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0580          hitsetitr != hitsetrange.second;
0581          ++hitsetitr)
0582     {
0583       // if(count>0)continue;
0584       TrkrHitSet* hitset = hitsetitr->second;
0585       unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0586       int side = TpcDefs::getSide(hitsetitr->first);
0587       unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0588 
0589       TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0590 
0591       for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0592            hitr != hitrangei.second;
0593            ++hitr)
0594       {
0595         unsigned short phibin = TpcDefs::getPad(hitr->first);
0596         unsigned short tbin = TpcDefs::getTBin(hitr->first);
0597         unsigned short adc = (hitr->second->getAdc());
0598 
0599         unsigned int pad_key = create_pad_key(side, layer, phibin);
0600 
0601         float fee = 0;
0602         std::map<unsigned int, chan_info>::iterator chan_it = chan_map.find(pad_key);
0603         if (chan_it != chan_map.end())
0604         {
0605           chan_info cinfo = (*chan_it).second;
0606           fee = cinfo.fee;
0607           // hpedestal2 = cinfo.ped;
0608           // hpedwidth2 = cinfo.width;
0609         }
0610 
0611         int rx = get_rx(layer);
0612         float corr = 0;
0613 
0614         unsigned int fee_key = create_fee_key(side, sector, rx, fee);
0615         std::map<unsigned int, std::vector<float>>::iterator fee_blm_it = feebaseline_map.find(fee_key);
0616         if (fee_blm_it != feebaseline_map.end())
0617         {
0618           if (tbin < (int) (*fee_blm_it).second.size())
0619           {
0620             corr = (*fee_blm_it).second[tbin];
0621           }
0622           hitr->second->setAdc(0);
0623           float nuadc = (float(adc) - corr);
0624           nuadc = std::max<float>(nuadc, 0);
0625           hitr->second->setAdc(nuadc);
0626 
0627           if (m_writeTree)
0628           {
0629             float fXh[18];
0630             int nh = 0;
0631 
0632             fXh[nh++] = _ievent - 1;
0633             fXh[nh++] = 0;       // gtm_bco;
0634             fXh[nh++] = 0;       // packet_id;
0635             fXh[nh++] = 0;       // ep;
0636             fXh[nh++] = sector;  // mc_sectors[sector % 12];//Sector;
0637             fXh[nh++] = side;
0638             fXh[nh++] = fee;
0639             fXh[nh++] = 0;  // channel;
0640             fXh[nh++] = 0;  // sampadd;
0641             fXh[nh++] = 0;  // sampch;
0642             fXh[nh++] = (float) phibin;
0643             fXh[nh++] = (float) tbin;
0644             fXh[nh++] = layer;
0645             fXh[nh++] = float(adc);
0646             fXh[nh++] = 0;  // hpedestal2;
0647             fXh[nh++] = 0;  // hpedwidth2;
0648             fXh[nh++] = corr;
0649 
0650             m_ntup_hits_corr->Fill(fXh);
0651           }
0652         }
0653       }
0654     }
0655   }
0656   // reset histogramms
0657   for (auto& hiter2 : feeadc_map)
0658   {
0659     if (hiter2.second != nullptr)
0660     {
0661       hiter2.second->Reset();
0662     }
0663   }
0664   feebaseline_map.clear();
0665   for (auto& hiter_entries : feeentries_map)
0666   {
0667     hiter_entries.second.assign(hiter_entries.second.size(), 0);
0668   }
0669 
0670   if (Verbosity())
0671   {
0672     std::cout << " event BCO: " << bco_min << " - " << bco_max << std::endl;
0673     std::cout << "TpcCombinedRawDataUnpacker:: done" << std::endl;
0674   }
0675 
0676   return Fun4AllReturnCodes::EVENT_OK;
0677 }
0678 
0679 int TpcCombinedRawDataUnpacker::End(PHCompositeNode* /*topNode*/)
0680 {
0681   if (m_writeTree)
0682   {
0683     m_file->cd();
0684     m_ntup->Write();
0685     m_ntup_hits->Write();
0686     m_ntup_hits_corr->Write();
0687     if (m_doChanHitsCut)
0688     {
0689       m_HitsinChan->Write();
0690     }
0691     m_file->Close();
0692   }
0693   if (Verbosity())
0694   {
0695     std::cout << "TpcCombinedRawDataUnpacker::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0696   }
0697   // if(m_Debug==1) hm->dumpHistos(m_filename, "RECREATE");
0698 
0699   return Fun4AllReturnCodes::EVENT_OK;
0700 }
[ 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!