Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:16:11 

0001 /*
0002  * TPCFEETestRecov1.cc
0003  *
0004  *  Created on: Sep 19, 2018
0005  *      Author: jinhuang
0006  */
0007 
0008 #include "TPCFEETestRecov1.h"
0009 
0010 #include "TPCDaqDefs.h"
0011 
0012 #include <g4detectors/PHG4Cell.h>
0013 #include <g4detectors/PHG4CellContainer.h>
0014 #include <g4detectors/PHG4CylinderCellGeom.h>
0015 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0016 #include <g4main/PHG4Hit.h>
0017 #include <g4main/PHG4HitContainer.h>
0018 
0019 #include <fun4all/Fun4AllHistoManager.h>
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021 #include <fun4all/Fun4AllServer.h>
0022 #include <fun4all/PHTFileServer.h>
0023 
0024 #include <phhepmc/PHHepMCGenEvent.h>
0025 #include <phhepmc/PHHepMCGenEventMap.h>
0026 #include <phool/PHCompositeNode.h>
0027 #include <phool/getClass.h>
0028 
0029 #include <Event/Event.h>
0030 #include <Event/EventTypes.h>
0031 #include <Event/packet.h>
0032 //#include <Event/packetConstants.h>
0033 #include <Event/oncsSubConstants.h>
0034 
0035 #include <TClonesArray.h>
0036 #include <TFile.h>
0037 #include <TH1D.h>
0038 #include <TH2D.h>
0039 #include <TString.h>
0040 #include <TTree.h>
0041 #include <TVector3.h>
0042 
0043 #include <CLHEP/Units/SystemOfUnits.h>
0044 
0045 #include <boost/bimap.hpp>
0046 #include <boost/bind.hpp>
0047 #include <boost/format.hpp>
0048 #include <boost/graph/adjacency_list.hpp>
0049 #include <boost/graph/connected_components.hpp>
0050 
0051 #include <algorithm>
0052 #include <array>
0053 #include <cassert>
0054 #include <cmath>
0055 #include <iostream>
0056 #include <limits>
0057 #include <map>
0058 #include <sstream>
0059 #include <stdexcept>
0060 #include <tuple>
0061 
0062 using namespace std;
0063 using namespace TPCDaqDefs::FEEv1;
0064 
0065 TPCFEETestRecov1::TPCFEETestRecov1(const std::string& outputfilename)
0066   : SubsysReco("TPCFEETestRecov1")
0067   , m_outputFileName(outputfilename)
0068   , m_eventT(nullptr)
0069   , m_peventHeader(&m_eventHeader)
0070   , m_nClusters(-1)
0071   , m_IOClusters(nullptr)
0072   , m_chanT(nullptr)
0073   , m_pchanHeader(&m_chanHeader)
0074   , m_chanData(kSAMPLE_LENGTH, 0)
0075   , m_clusteringZeroSuppression(50)
0076   , m_nPreSample(5)
0077   , m_nPostSample(5)
0078   , m_XRayLocationX(-1)
0079   , m_XRayLocationY(-1)
0080   , m_pdfMaker(nullptr)
0081 {
0082 }
0083 
0084 TPCFEETestRecov1::~TPCFEETestRecov1()
0085 {
0086   if (m_IOClusters)
0087   {
0088     m_IOClusters->Clear();
0089     delete m_IOClusters;
0090   }
0091 
0092   if (m_pdfMaker)
0093   {
0094     delete m_pdfMaker;
0095   }
0096 }
0097 
0098 int TPCFEETestRecov1::ResetEvent(PHCompositeNode* topNode)
0099 {
0100   m_eventHeader = EventHeader();
0101   m_padPlaneData.Reset();
0102   m_clusters.clear();
0103   m_chanHeader = ChannelHeader();
0104 
0105   m_nClusters = -1;
0106   assert(m_IOClusters);
0107   m_IOClusters->Clear();
0108 
0109   return Fun4AllReturnCodes::EVENT_OK;
0110 }
0111 
0112 int TPCFEETestRecov1::Init(PHCompositeNode* topNode)
0113 {
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116 
0117 int TPCFEETestRecov1::InitRun(PHCompositeNode* topNode)
0118 {
0119   if (Verbosity() >= VERBOSITY_SOME)
0120   {
0121     cout << "TPCFEETestRecov1::get_HistoManager - Making PHTFileServer " << m_outputFileName
0122          << endl;
0123 
0124     m_pdfMaker = new SampleFit_PowerLawDoubleExp_PDFMaker();
0125   }
0126   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0127 
0128   Fun4AllHistoManager* hm = getHistoManager();
0129   assert(hm);
0130 
0131   TH1D* h = new TH1D("hNormalization",  //
0132                      "Normalization;Items;Summed quantity", 10, .5, 10.5);
0133   int i = 1;
0134   h->GetXaxis()->SetBinLabel(i++, "Event count");
0135   h->GetXaxis()->SetBinLabel(i++, "Collision count");
0136   h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
0137   h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
0138   h->GetXaxis()->SetBinLabel(i++, "TPC Pad Hit");
0139   h->GetXaxis()->SetBinLabel(i++, "TPC Charge e");
0140   h->GetXaxis()->SetBinLabel(i++, "TPC Charge fC");
0141   h->GetXaxis()->LabelsOption("v");
0142   hm->registerHisto(h);
0143 
0144   m_eventT = new TTree("eventT", "TPC FEE per-event Tree");
0145   assert(m_eventT);
0146   m_eventT->Branch("evthdr", &m_peventHeader);
0147   m_eventT->Branch("nClusters", &m_nClusters, "nClusters/I");
0148   m_IOClusters = new TClonesArray("TPCFEETestRecov1::ClusterData", 1000);
0149   m_eventT->Branch("Clusters", &m_IOClusters);
0150 
0151   m_chanT = new TTree("chanT", "TPC FEE per-channel Tree");
0152   assert(m_chanT);
0153   m_chanT->Branch("event", &m_eventHeader.event, "event/I");
0154   m_chanT->Branch("chanhdr", &m_pchanHeader);
0155   m_chanT->Branch("adc", m_chanData.data(), str(boost::format("adc[%d]/i") % kSAMPLE_LENGTH).c_str());
0156 
0157   //  for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0158   //  {
0159   //    const PHG4CylinderCellGeom* layer_geom = seggeo->GetLayerCellGeom(layer);
0160 
0161   //    const string histNameCellHit(boost::str(boost::format{"hCellHit_Layer%1%"} % layer));
0162   //    const string histNameCellCharge(boost::str(boost::format{"hCellCharge_Layer%1%"} % layer));
0163 
0164   //  }
0165 
0166   //  hm->registerHisto(new TH2D("hLayerCellHit",  //
0167   //                             "Number of ADC time-bin hit per channel;Layer ID;Hit number",
0168   //                             m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0169   //                             300, -.5, 299.5));
0170   //  hm->registerHisto(new TH2D("hLayerCellCharge",  //
0171   //                             "Charge integrated over drift window per channel;Layer ID;Charge [fC]",
0172   //                             m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0173   //                             1000, 0, 1e7 * eplus / (1e-15 * coulomb)));
0174   //
0175   //  hm->registerHisto(new TH2D("hLayerSumCellHit",  //
0176   //                             "Number of ADC time-bin hit integrated over channels per layer;Layer ID;Hit number",
0177   //                             m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0178   //                             10000, -.5, 99999.5));
0179   //  hm->registerHisto(new TH2D("hLayerSumCellCharge",  //
0180   //                             "Charge integrated over drift window and channel per layer;Layer ID;Charge [fC]",
0181   //                             m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0182   //                             10000, 0, 1000 * 4e6 * eplus / (1e-15 * coulomb)));
0183 
0184   return Fun4AllReturnCodes::EVENT_OK;
0185 }
0186 
0187 int TPCFEETestRecov1::End(PHCompositeNode* topNode)
0188 {
0189   if (Verbosity() >= VERBOSITY_SOME)
0190   {
0191     cout << "TPCFEETestRecov1::End - write to " << m_outputFileName << endl;
0192   }
0193   PHTFileServer::get().cd(m_outputFileName);
0194 
0195   Fun4AllHistoManager* hm = getHistoManager();
0196   assert(hm);
0197   for (unsigned int i = 0; i < hm->nHistos(); i++)
0198     hm->getHisto(i)->Write();
0199 
0200   // help index files with TChain
0201   TTree* T_Index = new TTree("T_Index", "T_Index");
0202   assert(T_Index);
0203   T_Index->Write();
0204 
0205   m_eventT->Write();
0206   m_chanT->Write();
0207 
0208   if (m_pdfMaker)
0209   {
0210     delete m_pdfMaker;
0211     m_pdfMaker = nullptr;
0212   }
0213   return Fun4AllReturnCodes::EVENT_OK;
0214 }
0215 
0216 int TPCFEETestRecov1::process_event(PHCompositeNode* topNode)
0217 {
0218   Fun4AllHistoManager* hm = getHistoManager();
0219   assert(hm);
0220   TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
0221   assert(h_norm);
0222 
0223   Event* event = findNode::getClass<Event>(topNode, "PRDF");
0224   if (event == nullptr)
0225   {
0226     if (Verbosity() >= VERBOSITY_SOME)
0227       cout << "GenericUnpackPRDF::Process_Event - Event not found" << endl;
0228     return Fun4AllReturnCodes::DISCARDEVENT;
0229   }
0230 
0231   if (Verbosity() >= VERBOSITY_SOME)
0232     event->identify();
0233 
0234   // search for data event
0235   if (event->getEvtType() == BEGRUNEVENT)
0236   {
0237     get_motor_loc(event);
0238 
0239     return Fun4AllReturnCodes::EVENT_OK;
0240   }
0241   if (event->getEvtType() != DATAEVENT)
0242     return Fun4AllReturnCodes::DISCARDEVENT;
0243 
0244   m_eventHeader.run = event->getRunNumber();
0245   m_eventHeader.event = event->getEvtSequence();
0246 
0247   m_eventHeader.xray_x = m_XRayLocationX;
0248   m_eventHeader.xray_y = m_XRayLocationY;
0249 
0250   if (m_pdfMaker)
0251   {
0252     m_pdfMaker->MakeSectionPage(str(boost::format("ADC signal fits for Run %1% and event %2%") % m_eventHeader.run % m_eventHeader.event));
0253   }
0254 
0255   Packet* p = event->getPacket(kPACKET_ID, ID4EVT);
0256   if (p == nullptr)
0257     return Fun4AllReturnCodes::DISCARDEVENT;
0258 
0259   if (Verbosity() >= VERBOSITY_SOME) p->identify();
0260 
0261   if (Verbosity() >= VERBOSITY_MORE)
0262   {
0263     cout << "TPCFEETestRecov1::process_event - p->iValue(0) = "
0264          << p->iValue(0) << ", p->iValue(1) = " << p->iValue(1)
0265          << ", p->iValue(2) = " << p->iValue(2)
0266          << ", p->iValue(3) = " << p->iValue(3) << endl;
0267     p->dump();
0268   }
0269 
0270   m_eventHeader.bx_counter = 0;
0271   bool first_channel = true;
0272   for (unsigned int channel = 0; channel < kN_CHANNELS; channel++)
0273   {
0274     m_chanHeader = ChannelHeader();
0275 
0276     m_chanHeader.size = p->iValue(channel * kPACKET_LENGTH + 1) & 0xffff;         // number of words until the next channel (header included). this is the real packet_length
0277     m_chanHeader.packet_type = p->iValue(channel * kPACKET_LENGTH + 2) & 0xffff;  // that's the Elink packet type
0278     m_chanHeader.bx_counter = ((p->iValue(channel * kPACKET_LENGTH + 4) & 0xffff) << 4) | (p->iValue(channel * kPACKET_LENGTH + 5) & 0xffff);
0279     m_chanHeader.sampa_address = (p->iValue(channel * kPACKET_LENGTH + 3) >> 5) & 0xf;
0280     m_chanHeader.sampa_channel = p->iValue(channel * kPACKET_LENGTH + 3) & 0x1f;
0281     m_chanHeader.fee_channel = (m_chanHeader.sampa_address << 5) | m_chanHeader.sampa_channel;
0282 
0283     const pair<int, int> pad = SAMPAChan2PadXY(m_chanHeader.fee_channel);
0284 
0285     m_chanHeader.pad_x = pad.first;
0286     m_chanHeader.pad_y = pad.second;
0287 
0288     if (first_channel)
0289     {
0290       first_channel = false;
0291       m_eventHeader.bx_counter = m_chanHeader.bx_counter;
0292     }
0293     else if (m_eventHeader.bx_counter != m_chanHeader.bx_counter)
0294     {
0295       m_eventHeader.bx_counter_consistent = false;
0296 
0297       //      printf("TPCFEETestRecov1::process_event - ERROR: Malformed packet, event number %i, reason: bx_counter mismatch (expected 0x%x, got 0x%x)\n", m_eventHeader.event, m_eventHeader.bx_counter, m_chanHeader.bx_counter);
0298       //
0299       //      event->identify();
0300       //      p->identify();
0301       //      return Fun4AllReturnCodes::DISCARDEVENT;
0302     }
0303 
0304     if (m_chanHeader.fee_channel > 255 || m_chanHeader.sampa_address > 7 || m_chanHeader.sampa_channel > 31)
0305     {
0306       printf("TPCFEETestRecov1::process_event - ERROR: Malformed packet, event number %i, reason: bad channel (got %i, sampa_addr: %i, sampa_chan: %i)\n", m_eventHeader.event, m_chanHeader.fee_channel, m_chanHeader.sampa_address, m_chanHeader.sampa_channel);
0307 
0308       event->identify();
0309       p->identify();
0310       return Fun4AllReturnCodes::DISCARDEVENT;
0311     }
0312 
0313     //    SampaChannel *chan = fee_data->append(new SampaChannel(fee_channel, bx_counter, packet_type));
0314 
0315     assert(m_chanData.size() == kSAMPLE_LENGTH);
0316     fill(m_chanData.begin(), m_chanData.end(), 0);
0317     for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0318     {
0319       //        chan->append(p->iValue(channel * PACKET_LENGTH + 9 + sample) & 0xffff);
0320       uint32_t value = p->iValue(channel * kPACKET_LENGTH + 9 + sample) & 0xffff;
0321       m_chanData[sample] = value;
0322     }
0323 
0324     if (Verbosity() >= VERBOSITY_MORE)
0325     {
0326       cout << "TPCFEETestRecov1::process_event - "
0327            << "m_chanHeader.m_size = " << int(m_chanHeader.size) << ", "
0328            << "m_chanHeader.m_packet_type = " << int(m_chanHeader.packet_type) << ", "
0329            << "m_chanHeader.m_bx_counter = " << int(m_chanHeader.bx_counter) << ", "
0330            << "m_chanHeader.m_sampa_address = " << int(m_chanHeader.sampa_address) << ", "
0331            << "m_chanHeader.m_sampa_channel = " << int(m_chanHeader.sampa_channel) << ", "
0332            << "m_chanHeader.m_fee_channel = " << int(m_chanHeader.fee_channel) << ": "
0333            << " ";
0334 
0335       for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0336       {
0337         cout << "data[" << sample << "] = " << int(m_chanData[sample]) << " ";
0338       }
0339 
0340       cout << endl;
0341     }
0342 
0343     // fill event data
0344     if (PadPlaneData::IsValidPad(m_chanHeader.pad_x, m_chanHeader.pad_y))
0345     {
0346       vector<int>& paddata = m_padPlaneData.getPad(m_chanHeader.pad_x, m_chanHeader.pad_y);
0347 
0348       for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0349       {
0350         paddata[sample] = int(m_chanData[sample]);
0351       }
0352 
0353       auto pedestal_max = roughZeroSuppression(paddata);
0354       m_chanHeader.pedestal = pedestal_max.first;
0355       m_chanHeader.max = pedestal_max.second;
0356     }
0357     // output per-channel TTree
0358     m_chanT->Fill();
0359   }
0360 
0361   Clustering();
0362 
0363   h_norm->Fill("Event count", 1);
0364   m_eventT->Fill();
0365 
0366   return Fun4AllReturnCodes::EVENT_OK;
0367 }
0368 
0369 void TPCFEETestRecov1::Clustering()
0370 {
0371   // find cluster
0372   m_padPlaneData.Clustering(m_clusteringZeroSuppression, Verbosity() >= VERBOSITY_SOME);
0373   const multimap<int, PadPlaneData::SampleID>& groups = m_padPlaneData.getGroups();
0374 
0375   // export clusters
0376   assert(m_clusters.size() == 0);  //already cleared.
0377   for (const auto& iter : groups)
0378   {
0379     const int& i = iter.first;
0380     const PadPlaneData::SampleID& id = iter.second;
0381     m_clusters[i].padxs.insert(id.padx);
0382     m_clusters[i].padys.insert(id.pady);
0383     m_clusters[i].samples.insert(id.sample);
0384   }
0385 
0386   // process cluster
0387   for (auto& iter : m_clusters)
0388   {
0389     ClusterData& cluster = iter.second;
0390 
0391     assert(cluster.padxs.size() > 0);
0392     assert(cluster.padys.size() > 0);
0393     assert(cluster.samples.size() > 0);
0394 
0395     cluster.min_sample = max(0, *cluster.samples.begin() - m_nPreSample);
0396     cluster.max_sample = min((int) (kSAMPLE_LENGTH) -1, *cluster.samples.rbegin() + m_nPostSample);
0397     const int n_sample = cluster.max_sample - cluster.min_sample + 1;
0398 
0399     cluster.sum_samples.assign(n_sample, 0);
0400     for (int pad_x = *cluster.padxs.begin(); pad_x <= *cluster.padxs.rbegin(); ++pad_x)
0401     {
0402       cluster.padx_samples[pad_x].assign(n_sample, 0);
0403     }
0404     for (int pad_y = *cluster.padys.begin(); pad_y <= *cluster.padys.rbegin(); ++pad_y)
0405     {
0406       cluster.pady_samples[pad_y].assign(n_sample, 0);
0407     }
0408 
0409     for (int pad_x = *cluster.padxs.begin(); pad_x <= *cluster.padxs.rbegin(); ++pad_x)
0410     {
0411       for (int pad_y = *cluster.padys.begin(); pad_y <= *cluster.padys.rbegin(); ++pad_y)
0412       {
0413         assert(m_padPlaneData.IsValidPad(pad_x, pad_y));
0414 
0415         vector<int>& padsamples = m_padPlaneData.getPad(pad_x, pad_y);
0416 
0417         for (int i = 0; i < n_sample; ++i)
0418         {
0419           int adc = padsamples.at(cluster.min_sample + i);
0420           cluster.sum_samples[i] += adc;
0421           cluster.padx_samples[pad_x][i] += adc;
0422           cluster.pady_samples[pad_y][i] += adc;
0423         }
0424 
0425       }  //         for (int pad_y = *cluster.padys.begin(); pad_y<=*cluster.padys.rbegin() ;++pady)
0426 
0427     }  //    for (int pad_x = *cluster.padxs.begin(); pad_x<=*cluster.padxs.rbegin() ;++padx)
0428 
0429     if (m_pdfMaker)
0430     {
0431       m_pdfMaker->MakeSectionPage(str(boost::format("Event %1% Cluster %2%: sum all channel fit followed by fit of X/Y components") % m_eventHeader.event % iter.first));
0432     }
0433 
0434     // fit - overal cluster
0435     map<int, double> parameters_constraints;
0436     {
0437       double peak = NAN;
0438       double peak_sample = NAN;
0439       double pedstal = NAN;
0440       map<int, double> parameters_io;
0441       SampleFit_PowerLawDoubleExp(cluster.sum_samples, peak,
0442                                   peak_sample, pedstal, parameters_io, Verbosity());
0443 
0444       parameters_constraints[1] = parameters_io[1];
0445       parameters_constraints[2] = parameters_io[2];
0446       parameters_constraints[3] = parameters_io[3];
0447       parameters_constraints[5] = parameters_io[5];
0448       parameters_constraints[6] = parameters_io[6];
0449 
0450       cluster.peak = peak;
0451       cluster.peak_sample = peak_sample;
0452       cluster.pedstal = pedstal;
0453     }
0454 
0455     // fit - X
0456     {
0457       double sum_peak = 0;
0458       double sum_peak_padx = 0;
0459       for (int pad_x = *cluster.padxs.begin(); pad_x <= *cluster.padxs.rbegin(); ++pad_x)
0460       {
0461         double peak = NAN;
0462         double peak_sample = NAN;
0463         double pedstal = NAN;
0464         map<int, double> parameters_io(parameters_constraints);
0465 
0466         SampleFit_PowerLawDoubleExp(cluster.padx_samples[pad_x], peak,
0467                                     peak_sample, pedstal, parameters_io, Verbosity());
0468 
0469         cluster.padx_peaks[pad_x] = peak;
0470         sum_peak += peak;
0471         sum_peak_padx += peak * pad_x;
0472       }
0473       cluster.avg_padx = sum_peak_padx / sum_peak;
0474       cluster.size_pad_x = cluster.padxs.size();
0475     }
0476 
0477     // fit - Y
0478     {
0479       double sum_peak = 0;
0480       double sum_peak_pady = 0;
0481       for (int pad_y = *cluster.padys.begin(); pad_y <= *cluster.padys.rbegin(); ++pad_y)
0482       {
0483         double peak = NAN;
0484         double peak_sample = NAN;
0485         double pedstal = NAN;
0486         map<int, double> parameters_io(parameters_constraints);
0487 
0488         SampleFit_PowerLawDoubleExp(cluster.pady_samples[pad_y], peak,
0489                                     peak_sample, pedstal, parameters_io, Verbosity());
0490 
0491         cluster.pady_peaks[pad_y] = peak;
0492         sum_peak += peak;
0493         sum_peak_pady += peak * pad_y;
0494       }
0495       cluster.avg_pady = sum_peak_pady / sum_peak;
0496       cluster.size_pad_y = cluster.padys.size();
0497     }
0498   }  //   for (auto& iter : m_clusters)
0499 
0500   // sort by energy
0501   map<double, int> cluster_energy;
0502   for (auto& iter : m_clusters)
0503   {
0504     //reverse energy sorting
0505     cluster_energy[-iter.second.peak] = iter.first;
0506   }
0507 
0508   // save clusters
0509   m_nClusters = 0;
0510   assert(m_IOClusters);
0511   for (const auto& iter : cluster_energy)
0512   {
0513     ClusterData& cluster = m_clusters[iter.second];
0514 
0515     // super awkward ways of ROOT filling TClonesArray
0516     new ((*m_IOClusters)[m_nClusters++]) ClusterData(cluster);
0517   }
0518 }
0519 
0520 TPCFEETestRecov1::PadPlaneData::
0521     PadPlaneData()
0522   : m_data(kMaxPadY, vector<vector<int>>(kMaxPadX, vector<int>(kSAMPLE_LENGTH, 0)))
0523 {
0524 }
0525 
0526 void TPCFEETestRecov1::PadPlaneData::Reset()
0527 {
0528   for (auto& padrow : m_data)
0529   {
0530     for (auto& pad : padrow)
0531     {
0532       fill(pad.begin(), pad.end(), 0);
0533     }
0534   }
0535 
0536   m_groups.clear();
0537 }
0538 
0539 bool TPCFEETestRecov1::PadPlaneData::IsValidPad(const int pad_x, const int pad_y)
0540 {
0541   return (pad_x >= 0) and
0542          (pad_x < int(kMaxPadX)) and
0543          (pad_y >= 0) and
0544          (pad_y < int(kMaxPadY));
0545 }
0546 
0547 vector<int>& TPCFEETestRecov1::PadPlaneData::getPad(const int pad_x, const int pad_y)
0548 {
0549   assert(pad_x >= 0);
0550   assert(pad_x < int(kMaxPadX));
0551   assert(pad_y >= 0);
0552   assert(pad_y < int(kMaxPadY));
0553 
0554   return m_data[pad_y][pad_x];
0555 }
0556 
0557 std::pair<int, int> TPCFEETestRecov1::roughZeroSuppression(std::vector<int>& data)
0558 {
0559   std::vector<int> sorted_data(data);
0560 
0561   sort(sorted_data.begin(), sorted_data.end());
0562 
0563   const int pedestal = sorted_data[sorted_data.size() / 2];
0564   const int max = sorted_data.back();
0565 
0566   for (auto& d : data)
0567     d -= pedestal;
0568 
0569   return make_pair(pedestal, max);
0570 }
0571 
0572 bool operator<(const TPCFEETestRecov1::PadPlaneData::SampleID& s1, const TPCFEETestRecov1::PadPlaneData::SampleID& s2)
0573 {
0574   if (s1.pady == s2.pady)
0575   {
0576     if (s1.padx == s2.padx)
0577     {
0578       return s1.sample < s2.sample;
0579     }
0580     else
0581       return s1.padx < s2.padx;
0582   }
0583   else
0584     return s1.pady < s2.pady;
0585 }
0586 
0587 //! 3-D Graph clustering based on PHMakeGroups()
0588 void TPCFEETestRecov1::PadPlaneData::Clustering(int zero_suppression, bool verbosity)
0589 {
0590   using namespace boost;
0591   typedef adjacency_list<vecS, vecS, undirectedS> Graph;
0592   typedef bimap<Graph::vertex_descriptor, SampleID> VertexList;
0593 
0594   Graph G;
0595   VertexList vertex_list;
0596 
0597   for (unsigned int pady = 0; pady < kMaxPadY; ++pady)
0598   {
0599     for (unsigned int padx = 0; padx < kMaxPadX; ++padx)
0600     {
0601       for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0602       {
0603         if (m_data[pady][padx][sample] > zero_suppression)
0604         {
0605           SampleID id{(int) (pady), (int) (padx), (int) (sample)};
0606           Graph::vertex_descriptor v = boost::add_vertex(G);
0607           vertex_list.insert(VertexList::value_type(v, id));
0608 
0609           add_edge(v, v, G);
0610         }
0611       }  //      for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0612     }
0613   }  //   for (unsigned int pady = 0; pady < kMaxPadY; ++pady)
0614 
0615   // connect 3-D adjacent samples
0616   vector<SampleID> search_directions;
0617   search_directions.push_back(SampleID{0, 0, 1});
0618   search_directions.push_back(SampleID{0, 1, 0});
0619   search_directions.push_back(SampleID{1, 0, 0});
0620 
0621   for (const auto& it : vertex_list.right)
0622   {
0623     const SampleID id = it.first;
0624     const Graph::vertex_descriptor v = it.second;
0625 
0626     for (const SampleID& search_direction : search_directions)
0627     {
0628       //      const SampleID next_id = id + search_direction;
0629       SampleID next_id(id);
0630       next_id.adjust(search_direction);
0631 
0632       auto next_it = vertex_list.right.find(next_id);
0633       if (next_it != vertex_list.right.end())
0634       {
0635         add_edge(v, next_it->second, G);
0636       }
0637     }
0638 
0639   }  //  for (const auto & it : vertex_list)
0640 
0641   // Find the connections between the vertices of the graph (vertices are the rawhits,
0642   // connections are made when they are adjacent to one another)
0643   std::vector<int> component(num_vertices(G));
0644   connected_components(G, &component[0]);
0645 
0646   // Loop over the components(vertices) compiling a list of the unique
0647   // connections (ie clusters).
0648   set<int> comps;                // Number of unique components
0649   assert(m_groups.size() == 0);  // no overwrite
0650 
0651   for (unsigned int i = 0; i < component.size(); i++)
0652   {
0653     comps.insert(component[i]);
0654     m_groups.insert(make_pair(component[i], vertex_list.left.find(vertex(i, G))->second));
0655   }
0656 
0657   //debug prints
0658   if (verbosity)
0659     for (const int& comp : comps)
0660     {
0661       cout << "TPCFEETestRecov1::PadPlaneData::Clustering - find cluster " << comp << " containing ";
0662       const auto range = m_groups.equal_range(comp);
0663 
0664       for (auto iter = range.first; iter != range.second; ++iter)
0665       {
0666         const SampleID& id = iter->second;
0667         cout << "adc[" << id.pady << "][" << id.padx << "][" << id.sample << "] = " << m_data[id.pady][id.padx][id.sample] << ", ";
0668       }
0669       cout << endl;
0670     }  //  for (const int& comp : comps)
0671 }
0672 
0673 Fun4AllHistoManager*
0674 TPCFEETestRecov1::getHistoManager()
0675 {
0676   static string histname("TPCFEETestRecov1_HISTOS");
0677 
0678   Fun4AllServer* se = Fun4AllServer::instance();
0679   Fun4AllHistoManager* hm = se->getHistoManager(histname);
0680 
0681   if (not hm)
0682   {
0683     cout
0684         << "TPCFEETestRecov1::get_HistoManager - Making Fun4AllHistoManager "
0685         << histname << endl;
0686     hm = new Fun4AllHistoManager(histname);
0687     se->registerHistoManager(hm);
0688   }
0689 
0690   assert(hm);
0691 
0692   return hm;
0693 }
0694 
0695 void TPCFEETestRecov1::get_motor_loc(Event* evt)
0696 {
0697   assert(evt);
0698 
0699   Packet* motor_loc_p = evt->getPacket(910, IDCSTR);
0700 
0701   if (motor_loc_p)
0702   {
0703     string content;
0704 
0705     for (int i = 0; i < motor_loc_p->getLength(); i++)
0706     {
0707       content.push_back((char) motor_loc_p->iValue(i));
0708     }
0709 
0710     stringstream is(content);
0711     is >> m_XRayLocationX >> m_XRayLocationY;
0712 
0713     if (is.fail())
0714     {
0715       cout << "TPCFEETestRecov1::get_motor_loc - failed to load motor location from record [" << content << "]" << endl;
0716     }
0717     else if (Verbosity())
0718       cout << "TPCFEETestRecov1::get_motor_loc - received motor location " << m_XRayLocationX << ", " << m_XRayLocationY << " from record [" << content << "]" << endl;
0719   }
0720 }
[ 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!