29 #include <Event/Event.h>
30 #include <Event/EventTypes.h>
31 #include <Event/packet.h>
33 #include <Event/oncsSubConstants.h>
35 #include <TClonesArray.h>
43 #include <CLHEP/Units/SystemOfUnits.h>
45 #include <boost/bimap.hpp>
46 #include <boost/bind.hpp>
47 #include <boost/format.hpp>
48 #include <boost/graph/adjacency_list.hpp>
49 #include <boost/graph/connected_components.hpp>
63 using namespace TPCDaqDefs::FEEv1;
67 , m_outputFileName(outputfilename)
69 , m_peventHeader(&m_eventHeader)
71 , m_IOClusters(nullptr)
73 , m_pchanHeader(&m_chanHeader)
75 , m_clusteringZeroSuppression(50)
121 cout <<
"TPCFEETestRecov1::get_HistoManager - Making PHTFileServer " <<
m_outputFileName
131 TH1D*
h =
new TH1D(
"hNormalization",
132 "Normalization;Items;Summed quantity", 10, .5, 10.5);
134 h->GetXaxis()->SetBinLabel(i++,
"Event count");
135 h->GetXaxis()->SetBinLabel(i++,
"Collision count");
136 h->GetXaxis()->SetBinLabel(i++,
"TPC G4Hit");
137 h->GetXaxis()->SetBinLabel(i++,
"TPC G4Hit Edep");
138 h->GetXaxis()->SetBinLabel(i++,
"TPC Pad Hit");
139 h->GetXaxis()->SetBinLabel(i++,
"TPC Charge e");
140 h->GetXaxis()->SetBinLabel(i++,
"TPC Charge fC");
141 h->GetXaxis()->LabelsOption(
"v");
144 m_eventT =
new TTree(
"eventT",
"TPC FEE per-event Tree");
148 m_IOClusters =
new TClonesArray(
"TPCFEETestRecov1::ClusterData", 1000);
151 m_chanT =
new TTree(
"chanT",
"TPC FEE per-channel Tree");
197 for (
unsigned int i = 0;
i < hm->
nHistos();
i++)
201 TTree* T_Index =
new TTree(
"T_Index",
"T_Index");
220 TH1D* h_norm =
dynamic_cast<TH1D*
>(hm->
getHisto(
"hNormalization"));
223 Event*
event = findNode::getClass<Event>(topNode,
"PRDF");
224 if (
event ==
nullptr)
227 cout <<
"GenericUnpackPRDF::Process_Event - Event not found" << endl;
263 cout <<
"TPCFEETestRecov1::process_event - p->iValue(0) = "
265 <<
", p->iValue(2) = " << p->
iValue(2)
266 <<
", p->iValue(3) = " << p->
iValue(3) << endl;
271 bool first_channel =
true;
290 first_channel =
false;
326 cout <<
"TPCFEETestRecov1::process_event - "
337 cout <<
"data[" << sample <<
"] = " << int(
m_chanData[sample]) <<
" ";
363 h_norm->Fill(
"Event count", 1);
377 for (
const auto& iter : groups)
379 const int&
i = iter.first;
400 for (
int pad_x = *cluster.
padxs.begin(); pad_x <= *cluster.
padxs.rbegin(); ++pad_x)
404 for (
int pad_y = *cluster.
padys.begin(); pad_y <= *cluster.
padys.rbegin(); ++pad_y)
409 for (
int pad_x = *cluster.
padxs.begin(); pad_x <= *cluster.
padxs.rbegin(); ++pad_x)
411 for (
int pad_y = *cluster.
padys.begin(); pad_y <= *cluster.
padys.rbegin(); ++pad_y)
417 for (
int i = 0;
i < n_sample; ++
i)
435 map<int, double> parameters_constraints;
438 double peak_sample = NAN;
439 double pedstal = NAN;
440 map<int, double> parameters_io;
442 peak_sample, pedstal, parameters_io,
Verbosity());
444 parameters_constraints[1] = parameters_io[1];
445 parameters_constraints[2] = parameters_io[2];
446 parameters_constraints[3] = parameters_io[3];
447 parameters_constraints[5] = parameters_io[5];
448 parameters_constraints[6] = parameters_io[6];
458 double sum_peak_padx = 0;
459 for (
int pad_x = *cluster.
padxs.begin(); pad_x <= *cluster.
padxs.rbegin(); ++pad_x)
462 double peak_sample = NAN;
463 double pedstal = NAN;
464 map<int, double> parameters_io(parameters_constraints);
467 peak_sample, pedstal, parameters_io,
Verbosity());
471 sum_peak_padx += peak * pad_x;
473 cluster.
avg_padx = sum_peak_padx / sum_peak;
480 double sum_peak_pady = 0;
481 for (
int pad_y = *cluster.
padys.begin(); pad_y <= *cluster.
padys.rbegin(); ++pad_y)
484 double peak_sample = NAN;
485 double pedstal = NAN;
486 map<int, double> parameters_io(parameters_constraints);
489 peak_sample, pedstal, parameters_io,
Verbosity());
493 sum_peak_pady += peak * pad_y;
495 cluster.
avg_pady = sum_peak_pady / sum_peak;
501 map<double, int> cluster_energy;
502 for (
auto& iter : m_clusters)
505 cluster_energy[-iter.second.peak] = iter.first;
511 for (
const auto& iter : cluster_energy)
528 for (
auto& padrow : m_data)
530 for (
auto& pad : padrow)
532 fill(pad.begin(), pad.end(), 0);
541 return (pad_x >= 0) and
554 return m_data[pad_y][pad_x];
559 std::vector<int> sorted_data(data);
561 sort(sorted_data.begin(), sorted_data.end());
563 const int pedestal = sorted_data[sorted_data.size() / 2];
564 const int max = sorted_data.back();
569 return make_pair(pedestal, max);
590 using namespace boost;
591 typedef adjacency_list<vecS, vecS, undirectedS> Graph;
592 typedef bimap<Graph::vertex_descriptor, SampleID> VertexList;
595 VertexList vertex_list;
597 for (
unsigned int pady = 0; pady <
kMaxPadY; ++pady)
599 for (
unsigned int padx = 0; padx <
kMaxPadX; ++padx)
603 if (m_data[pady][padx][sample] > zero_suppression)
605 SampleID id{(int) (pady), (int) (padx), (int) (sample)};
606 Graph::vertex_descriptor
v = boost::add_vertex(
G);
616 vector<SampleID> search_directions;
617 search_directions.push_back(
SampleID{0, 0, 1});
618 search_directions.push_back(
SampleID{0, 1, 0});
619 search_directions.push_back(
SampleID{1, 0, 0});
621 for (
const auto&
it : vertex_list.right)
624 const Graph::vertex_descriptor
v =
it.second;
626 for (
const SampleID& search_direction : search_directions)
630 next_id.
adjust(search_direction);
632 auto next_it = vertex_list.right.find(next_id);
633 if (next_it != vertex_list.right.end())
635 add_edge(v, next_it->second,
G);
643 std::vector<int> component(num_vertices(
G));
644 connected_components(
G, &component[0]);
649 assert(m_groups.size() == 0);
651 for (
unsigned int i = 0;
i < component.size();
i++)
653 comps.insert(component[
i]);
654 m_groups.insert(make_pair(component[i], vertex_list.left.find(
vertex(i,
G))->second));
659 for (
const int&
comp : comps)
661 cout <<
"TPCFEETestRecov1::PadPlaneData::Clustering - find cluster " <<
comp <<
" containing ";
662 const auto range = m_groups.equal_range(comp);
664 for (
auto iter = range.first; iter != range.second; ++iter)
667 cout <<
"adc[" <<
id.
pady <<
"][" <<
id.padx <<
"][" <<
id.sample <<
"] = " << m_data[
id.pady][
id.padx][
id.sample] <<
", ";
676 static string histname(
"TPCFEETestRecov1_HISTOS");
684 <<
"TPCFEETestRecov1::get_HistoManager - Making Fun4AllHistoManager "
707 content.push_back((
char) motor_loc_p->
iValue(
i));
710 stringstream is(content);
715 cout <<
"TPCFEETestRecov1::get_motor_loc - failed to load motor location from record [" << content <<
"]" << endl;
718 cout <<
"TPCFEETestRecov1::get_motor_loc - received motor location " <<
m_XRayLocationX <<
", " << m_XRayLocationY <<
" from record [" << content <<
"]" << endl;