Back to home page

sPhenix code displayed by LXR
 
 

    

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

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/PHG4TpcGeom.h>
0012 #include <g4detectors/PHG4TpcGeomContainer.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   PHG4TpcGeomContainer* geom_container =
0246       findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0247   if (!geom_container)
0248   {
0249     std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << 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      
0309     int region = 2;
0310     if(layer < 7 + 16)
0311     {
0312       region = 0;
0313     }
0314     else if( layer < 7 + 32)
0315     {
0316       region = 1;
0317     }
0318 
0319     varname = "phi";  // + std::to_string(key);
0320     double phi = ((side == 1 ? 1 : -1) * (m_cdbttree->GetDoubleValue(key, varname) - M_PI / 2.)) + ((sector % 12) * M_PI / 6);
0321     PHG4TpcGeom* layergeom = geom_container->GetLayerCellGeom(layer);
0322     unsigned int phibin = layergeom->get_phibin(phi, side);
0323   
0324     hit_set_key = TpcDefs::genHitSetKey(layer, (mc_sectors[sector % 12]), side);
0325     hit_set_container_itr = trkr_hit_set_container->findOrAddHitSet(hit_set_key);
0326 
0327     float hpedestal = 0;
0328     float hpedwidth = 0;
0329 
0330     if (Verbosity() > 2)
0331     {
0332       std::cout << "TpcCombinedRawDataUnpacker:: do zero suppression" << std::endl;
0333     }
0334     TH2* feehist = nullptr;
0335     hpedestal = 60;
0336     hpedwidth = m_zs_threshold[region];
0337 
0338     unsigned int pad_key = create_pad_key(side, layer, phibin);
0339 
0340     std::map<unsigned int, chan_info>::iterator chan_it = chan_map.find(pad_key);
0341     if (chan_it != chan_map.end())
0342     {
0343       (*chan_it).second.ped = hpedestal;
0344       (*chan_it).second.width = hpedwidth;
0345     }
0346     else
0347     {
0348       chan_info nucinfo;
0349       nucinfo.fee = fee;
0350       nucinfo.ped = hpedestal;
0351       nucinfo.width = hpedwidth;
0352       chan_map.insert(std::make_pair(pad_key, nucinfo));
0353     }
0354     int rx = get_rx(layer);
0355     unsigned int fee_key = create_fee_key(side, mc_sectors[sector % 12], rx, fee);
0356     // find or insert TH2C;
0357     std::map<unsigned int, TH2*>::iterator fee_map_it;
0358 
0359     fee_map_it = feeadc_map.find(fee_key);
0360     if (fee_map_it != feeadc_map.end())
0361     {
0362       feehist = (*fee_map_it).second;
0363     }
0364     else
0365     {
0366       std::string histname = "h" + std::to_string(fee_key);
0367       feehist = new TH2C(histname.c_str(), "histname", max_time_range + 1, -0.5, max_time_range + 0.5, 501, -0.5, 1000.5);
0368       feeadc_map.insert(std::make_pair(fee_key, feehist));
0369       std::vector<int> feeentries(feehist->GetNbinsX(), 0);
0370       feeentries_map.insert(std::make_pair(fee_key, feeentries));
0371     }
0372     auto fee_entries_it = feeentries_map.find(fee_key);
0373     std::vector<int>& fee_entries_vec = (*fee_entries_it).second;
0374 
0375     float threshold_cut = m_zs_threshold[region];
0376 
0377     int nhitschan = 0;
0378 
0379     if (m_doChanHitsCut)
0380     {
0381       for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(tpchit->CreateAdcIterator()); !adc_iterator->IsDone(); adc_iterator->Next())
0382       {
0383         const uint16_t s = adc_iterator->CurrentTimeBin();
0384         const uint16_t adc = adc_iterator->CurrentAdc();
0385         int t = s - m_presampleShift - m_t0;
0386         if (t < 0)
0387         {
0388           continue;
0389         }
0390         if (feehist != nullptr)
0391         {
0392           if (adc > 0)
0393           {
0394             if ((float(adc) - hpedestal) > threshold_cut)
0395             {
0396               nhitschan++;
0397             }
0398           }
0399         }
0400       }
0401       if (m_writeTree)
0402       {
0403         m_HitChanDis->Fill(nhitschan, channel);
0404         m_HitsinChan->Fill(nhitschan);
0405       }
0406       if (nhitschan > m_ChanHitsCut)
0407       {
0408         continue;
0409       }
0410     }
0411 
0412     for (std::unique_ptr<TpcRawHit::AdcIterator> adc_iterator(tpchit->CreateAdcIterator());
0413          !adc_iterator->IsDone();
0414          adc_iterator->Next())
0415     {
0416       const uint16_t s = adc_iterator->CurrentTimeBin();
0417       const uint16_t adc = adc_iterator->CurrentAdc();
0418       int t = s - m_presampleShift - m_t0;
0419       if (t < 0)
0420       {
0421         continue;
0422       }
0423       if (feehist != nullptr)
0424       {
0425         if (adc > 0)
0426         {
0427           if ((float(adc) - hpedestal) > threshold_cut)
0428           {
0429             feehist->Fill(t, adc - hpedestal);
0430             if (t < (int) fee_entries_vec.size())
0431             {
0432               fee_entries_vec[t]++;
0433             }
0434           }
0435         }
0436       }
0437 
0438       if ((float(adc) - hpedestal) > threshold_cut)
0439       {
0440         hit_key = TpcDefs::genHitKey(phibin, (unsigned int) t);
0441         // find existing hit, or create new one
0442         hit = hit_set_container_itr->second->getHit(hit_key);
0443         if (!hit)
0444         {
0445           hit = new TrkrHitv2();
0446           hit->setAdc(float(adc) - hpedestal);
0447           hit_set_container_itr->second->addHitSpecificKey(hit_key, hit);
0448         }
0449 
0450         if (m_writeTree)
0451         {
0452           float fXh[18];
0453           int nh = 0;
0454 
0455           fXh[nh++] = _ievent - 1;
0456           fXh[nh++] = gtm_bco;                  // gtm_bco;
0457           fXh[nh++] = packet_id;                // packet_id;
0458           fXh[nh++] = ep;                       // ep;
0459           fXh[nh++] = mc_sectors[sector % 12];  // Sector;
0460           fXh[nh++] = side;
0461           fXh[nh++] = fee;
0462           fXh[nh++] = channel;  // channel;
0463           fXh[nh++] = sampadd;  // sampadd;
0464           fXh[nh++] = sampch;   // sampch;
0465           fXh[nh++] = (float) phibin;
0466           fXh[nh++] = (float) t;
0467           fXh[nh++] = layer;
0468           fXh[nh++] = (float(adc) - hpedestal);
0469           fXh[nh++] = hpedestal;
0470           fXh[nh++] = hpedwidth;
0471           m_ntup_hits->Fill(fXh);
0472         }
0473       }
0474     }
0475   }
0476 
0477   if (m_do_baseline_corr == true)
0478   {
0479     // Histos filled now process them for fee local baselines
0480 
0481     int nhistfilled = 0;
0482     int nhisttotal = 0;
0483     for (auto& hiter : feeadc_map)
0484     {
0485       if (hiter.second != nullptr)
0486       {
0487         unsigned int fee_key = hiter.first;
0488         unsigned int side;
0489         unsigned int sector;
0490         unsigned int rx;
0491         unsigned int fee;
0492         unpack_fee_key(side, sector, rx, fee, fee_key);
0493         TH2* hist2d = hiter.second;
0494         std::map<unsigned int, std::vector<int>>::iterator fee_entries_it = feeentries_map.find(fee_key);
0495         if (fee_entries_it == feeentries_map.end())
0496         {
0497           continue;
0498           //      fee_entries_vec_it = (*fee_entries_it).second.begin();
0499         }
0500         std::vector<int>::iterator fee_entries_vec_it = (*fee_entries_it).second.begin();
0501 
0502         std::vector<float> pedvec(hist2d->GetNbinsX(), 0);
0503         feebaseline_map.insert(std::make_pair(hiter.first, pedvec));
0504         std::map<unsigned int, std::vector<float>>::iterator fee_blm_it = feebaseline_map.find(hiter.first);
0505         (*fee_blm_it).second.resize(hist2d->GetNbinsX(), 0);
0506         for (int binx = 1; binx < hist2d->GetNbinsX(); binx++)
0507         {
0508           double timebin = (hist2d->GetXaxis())->GetBinCenter(binx);
0509           std::string histname1d = "h" + std::to_string(hiter.first) + "_" + std::to_string((int) timebin);
0510           nhisttotal++;
0511           float local_ped = 0;
0512           float local_width = 0;
0513           float entries = fee_entries_vec_it[timebin];
0514           if (fee_entries_vec_it[timebin] > 100)
0515           {
0516             nhistfilled++;
0517 
0518             TH1D* hist1d = hist2d->ProjectionY(histname1d.c_str(), binx, binx);
0519             if (hist1d->GetEntries() != fee_entries_vec_it[timebin])
0520             {
0521               std::cout << " vec " << fee_entries_vec_it[timebin]
0522                         << " hist " << hist1d->GetEntries()
0523                         << std::endl;
0524             }
0525             if (hist1d->GetEntries() > 10)
0526             {
0527               int maxbin = hist1d->GetMaximumBin();
0528               // calc peak position
0529               double hadc_sum = 0.0;
0530               double hibin_sum = 0.0;
0531               double hibin2_sum = 0.0;
0532 
0533               for (int isum = -3; isum <= 3; isum++)
0534               {
0535                 float val = hist1d->GetBinContent(maxbin + isum);
0536                 float center = hist1d->GetBinCenter(maxbin + isum);
0537                 hibin_sum += center * val;
0538                 hibin2_sum += center * center * val;
0539                 hadc_sum += val;
0540               }
0541               local_ped = hibin_sum / hadc_sum;
0542               local_width = sqrt((hibin2_sum / hadc_sum) - (local_ped * local_ped));
0543             }
0544             delete hist1d;
0545           }
0546           (*fee_blm_it).second[(int) timebin] = local_ped + m_baseline_nsigma * local_width;
0547 
0548           if (m_writeTree)
0549           {
0550             float fXh[11];
0551             int nh = 0;
0552 
0553             fXh[nh++] = _ievent - 1;
0554             fXh[nh++] = 0;                        // gtm_bco;
0555             fXh[nh++] = 0;                        // packet_id;
0556             fXh[nh++] = 0;                        // ep;
0557             fXh[nh++] = mc_sectors[sector % 12];  // Sector;
0558             fXh[nh++] = side;
0559             fXh[nh++] = fee;
0560             fXh[nh++] = rx;
0561             fXh[nh++] = entries;
0562             fXh[nh++] = local_ped;
0563             fXh[nh++] = local_width;
0564             m_ntup->Fill(fXh);
0565           }
0566         }
0567       }
0568     }
0569     if (Verbosity() >= 1)
0570     {
0571       std::cout << " filled " << nhistfilled
0572                 << " total " << nhisttotal
0573                 << std::endl;
0574 
0575       std::cout << "second loop " << m_do_baseline_corr << std::endl;
0576     }
0577 
0578     // second loop over hits to apply baseline correction
0579     TrkrHitSetContainer::ConstRange hitsetrange;
0580     hitsetrange = trkr_hit_set_container->getHitSets(TrkrDefs::TrkrId::tpcId);
0581 
0582     for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0583          hitsetitr != hitsetrange.second;
0584          ++hitsetitr)
0585     {
0586       // if(count>0)continue;
0587       TrkrHitSet* hitset = hitsetitr->second;
0588       unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0589       int side = TpcDefs::getSide(hitsetitr->first);
0590       unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0591 
0592       TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0593 
0594       for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0595            hitr != hitrangei.second;
0596            ++hitr)
0597       {
0598         unsigned short phibin = TpcDefs::getPad(hitr->first);
0599         unsigned short tbin = TpcDefs::getTBin(hitr->first);
0600         unsigned short adc = (hitr->second->getAdc());
0601 
0602         unsigned int pad_key = create_pad_key(side, layer, phibin);
0603 
0604         float fee = 0;
0605         std::map<unsigned int, chan_info>::iterator chan_it = chan_map.find(pad_key);
0606         if (chan_it != chan_map.end())
0607         {
0608           chan_info cinfo = (*chan_it).second;
0609           fee = cinfo.fee;
0610           // hpedestal2 = cinfo.ped;
0611           // hpedwidth2 = cinfo.width;
0612         }
0613 
0614         int rx = get_rx(layer);
0615         float corr = 0;
0616 
0617         unsigned int fee_key = create_fee_key(side, sector, rx, fee);
0618         std::map<unsigned int, std::vector<float>>::iterator fee_blm_it = feebaseline_map.find(fee_key);
0619         if (fee_blm_it != feebaseline_map.end())
0620         {
0621           if (tbin < (int) (*fee_blm_it).second.size())
0622           {
0623             corr = (*fee_blm_it).second[tbin];
0624           }
0625           hitr->second->setAdc(0);
0626           float nuadc = (float(adc) - corr);
0627           nuadc = std::max<float>(nuadc, 0);
0628           hitr->second->setAdc(nuadc);
0629 
0630           if (m_writeTree)
0631           {
0632             float fXh[18];
0633             int nh = 0;
0634 
0635             fXh[nh++] = _ievent - 1;
0636             fXh[nh++] = 0;       // gtm_bco;
0637             fXh[nh++] = 0;       // packet_id;
0638             fXh[nh++] = 0;       // ep;
0639             fXh[nh++] = sector;  // mc_sectors[sector % 12];//Sector;
0640             fXh[nh++] = side;
0641             fXh[nh++] = fee;
0642             fXh[nh++] = 0;  // channel;
0643             fXh[nh++] = 0;  // sampadd;
0644             fXh[nh++] = 0;  // sampch;
0645             fXh[nh++] = (float) phibin;
0646             fXh[nh++] = (float) tbin;
0647             fXh[nh++] = layer;
0648             fXh[nh++] = float(adc);
0649             fXh[nh++] = 0;  // hpedestal2;
0650             fXh[nh++] = 0;  // hpedwidth2;
0651             fXh[nh++] = corr;
0652 
0653             m_ntup_hits_corr->Fill(fXh);
0654           }
0655         }
0656       }
0657     }
0658   }
0659   // reset histogramms
0660   for (auto& hiter2 : feeadc_map)
0661   {
0662     if (hiter2.second != nullptr)
0663     {
0664       hiter2.second->Reset();
0665     }
0666   }
0667   feebaseline_map.clear();
0668   for (auto& hiter_entries : feeentries_map)
0669   {
0670     hiter_entries.second.assign(hiter_entries.second.size(), 0);
0671   }
0672 
0673   if (Verbosity())
0674   {
0675     std::cout << " event BCO: " << bco_min << " - " << bco_max << std::endl;
0676     std::cout << "TpcCombinedRawDataUnpacker:: done" << std::endl;
0677   }
0678 
0679   return Fun4AllReturnCodes::EVENT_OK;
0680 }
0681 
0682 int TpcCombinedRawDataUnpacker::End(PHCompositeNode* /*topNode*/)
0683 {
0684   if (m_writeTree)
0685   {
0686     m_file->cd();
0687     m_ntup->Write();
0688     m_ntup_hits->Write();
0689     m_ntup_hits_corr->Write();
0690     if (m_doChanHitsCut)
0691     {
0692       m_HitsinChan->Write();
0693     }
0694     m_file->Close();
0695   }
0696   if (Verbosity())
0697   {
0698     std::cout << "TpcCombinedRawDataUnpacker::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0699   }
0700   // if(m_Debug==1) hm->dumpHistos(m_filename, "RECREATE");
0701 
0702   return Fun4AllReturnCodes::EVENT_OK;
0703 }
[ 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!