packages/tpc/TpcCombinedRawDataUnpackerDebug.cc

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