2 #include "messagefacility/MessageLogger/MessageLogger.h"
24 mf::LogInfo(
"HitAnalysisAlg") <<
"HitAnalysisAlg configured\n";
41 fLocalDirName = pset.get<std::string>(
"LocalDirName", std::string(
"wow"));
47 TDirectory* rootDirectory)
60 fDriftTimes[0] = std::make_unique<TH1D>(
"DriftTime0",
";time(ticks)", 3200, 0., 9600.);
61 fDriftTimes[1] = std::make_unique<TH1D>(
"DriftTime1",
";time(ticks)", 3200, 0., 9600.);
62 fDriftTimes[2] = std::make_unique<TH1D>(
"DriftTime2",
";time(ticks)", 3200, 0., 9600.);
64 fHitsByTime[0] = std::make_unique<TH1D>(
"HitsByTime0",
";Tick", 1600, 0., 6400.);
65 fHitsByTime[1] = std::make_unique<TH1D>(
"HitsByTime1",
";Tick", 1600, 0., 6400.);
66 fHitsByTime[2] = std::make_unique<TH1D>(
"HitsByTime2",
";Tick", 1600, 0., 6400.);
68 fPulseHeight[0] = std::make_unique<TH1D>(
"PulseHeight0",
"PH (ADC)", 300, 0., 150.);
69 fPulseHeight[1] = std::make_unique<TH1D>(
"PulseHeight1",
"PH (ADC)", 300, 0., 150.);
70 fPulseHeight[2] = std::make_unique<TH1D>(
"PulseHeight2",
"PH (ADC)", 300, 0., 150.);
71 fPulseHeightSingle[0] = std::make_unique<TH1D>(
"PulseHeightS0",
"PH (ADC)", 300, 0., 150.);
72 fPulseHeightSingle[1] = std::make_unique<TH1D>(
"PulseHeightS1",
"PH (ADC)", 300, 0., 150.);
73 fPulseHeightSingle[2] = std::make_unique<TH1D>(
"PulseHeightS2",
"PH (ADC)", 300, 0., 150.);
74 fPulseHeightMulti[0] = std::make_unique<TH1D>(
"PulseHeightM0",
"PH (ADC)", 300, 0., 150.);
75 fPulseHeightMulti[1] = std::make_unique<TH1D>(
"PulseHeightM1",
"PH (ADC)", 300, 0., 150.);
76 fPulseHeightMulti[2] = std::make_unique<TH1D>(
"PulseHeightM2",
"PH (ADC)", 300, 0., 150.);
77 fChi2DOF[0] = std::make_unique<TH1D>(
"Chi2DOF0",
"Chi2DOF", 502, -1., 250.);
78 fChi2DOF[1] = std::make_unique<TH1D>(
"Chi2DOF1",
"Chi2DOF", 502, -1., 250.);
79 fChi2DOF[2] = std::make_unique<TH1D>(
"Chi2DOF2",
"Chi2DOF", 502, -1., 250.);
80 fNumDegFree[0] = std::make_unique<TH1D>(
"NumDegFree0",
"NDF", 100, 0., 100.);
81 fNumDegFree[1] = std::make_unique<TH1D>(
"NumDegFree1",
"NDF", 100, 0., 100.);
82 fNumDegFree[2] = std::make_unique<TH1D>(
"NumDegFree2",
"NDF", 100, 0., 100.);
83 fChi2DOFSingle[0] = std::make_unique<TH1D>(
"Chi2DOFS0",
"Chi2DOF", 502, -1., 250.);
84 fChi2DOFSingle[1] = std::make_unique<TH1D>(
"Chi2DOFS1",
"Chi2DOF", 502, -1., 250.);
85 fChi2DOFSingle[2] = std::make_unique<TH1D>(
"Chi2DOFS2",
"Chi2DOF", 502, -1., 250.);
86 fHitMult[0] = std::make_unique<TH1D>(
"HitMult0",
"# hits", 15, 0., 15.);
87 fHitMult[1] = std::make_unique<TH1D>(
"HitMult1",
"# hits", 15, 0., 15.);
88 fHitMult[2] = std::make_unique<TH1D>(
"HitMult2",
"# hits", 15, 0., 15.);
89 fHitCharge[0] = std::make_unique<TH1D>(
"HitCharge0",
"Charge", 1000, 0., 2000.);
90 fHitCharge[1] = std::make_unique<TH1D>(
"HitCharge1",
"Charge", 1000, 0., 2000.);
91 fHitCharge[2] = std::make_unique<TH1D>(
"HitCharge2",
"Charge", 1000, 0., 2000.);
92 fFitWidth[0] = std::make_unique<TH1D>(
"FitWidth0",
"Width", 100, 0., 10.);
93 fFitWidth[1] = std::make_unique<TH1D>(
"FitWidth1",
"Width", 100, 0., 10.);
94 fFitWidth[2] = std::make_unique<TH1D>(
"FitWidth2",
"Width", 100, 0., 10.);
95 fHitSumADC[0] = std::make_unique<TH1D>(
"SumADC0",
"Sum ADC", 1000, 0., 2000.);
96 fHitSumADC[1] = std::make_unique<TH1D>(
"SumADC1",
"Sum ADC", 1000, 0., 2000.);
97 fHitSumADC[2] = std::make_unique<TH1D>(
"SumADC2",
"Sum ADC", 1000, 0., 2000.);
99 fNDFVsChi2[0] = std::make_unique<TH2D>(
"NDFVsChi20",
";NDF;Chi2", 50, 0., 50., 101, -1., 100.);
100 fNDFVsChi2[1] = std::make_unique<TH2D>(
"NDFVsChi21",
";NDF;Chi2", 50, 0., 50., 101, -1., 100.);
101 fNDFVsChi2[2] = std::make_unique<TH2D>(
"NDFVsChi22",
";NDF;Chi2", 50, 0., 50., 101, -1., 100.);
103 fPulseHVsWidth[0] = std::make_unique<TH2D>(
"PHVsWidth0",
";PH;Width", 100, 0., 100., 100, 0., 20.);
104 fPulseHVsWidth[1] = std::make_unique<TH2D>(
"PHVsWidth1",
";PH;Width", 100, 0., 100., 100, 0., 20.);
105 fPulseHVsWidth[2] = std::make_unique<TH2D>(
"PHVsWidth2",
";PH;Width", 100, 0., 100., 100, 0., 20.);
107 fPulseHVsCharge[0] = std::make_unique<TH2D>(
"PHVsChrg0",
";PH;Q", 100, 0., 100., 100, 0., 2000.);
108 fPulseHVsCharge[1] = std::make_unique<TH2D>(
"PHVsChrg1",
";PH;Q", 100, 0., 100., 100, 0., 2000.);
109 fPulseHVsCharge[2] = std::make_unique<TH2D>(
"PHVsChrg2",
";PH;Q", 100, 0., 100., 100, 0., 2000.);
111 fPulseHVsHitNo[0] = std::make_unique<TProfile>(
"PHVsNo0",
";Hit #;PH", 1000, 0., 1000., 0., 100.);
112 fPulseHVsHitNo[1] = std::make_unique<TProfile>(
"PHVsNo1",
";Hit #;PH", 1000, 0., 1000., 0., 100.);
113 fPulseHVsHitNo[2] = std::make_unique<TProfile>(
"PHVsNo2",
";Hit #;PH", 1000, 0., 1000., 0., 100.);
115 fChargeVsHitNo[0] = std::make_unique<TProfile>(
"QVsNo0",
";Hit No;Q", 1000, 0., 1000., 0., 2000.);
116 fChargeVsHitNo[1] = std::make_unique<TProfile>(
"QVsNo1",
";Hit No;Q", 1000, 0., 1000., 0., 2000.);
117 fChargeVsHitNo[2] = std::make_unique<TProfile>(
"QVsNo2",
";Hit No;Q", 1000, 0., 1000., 0., 2000.);
119 fChargeVsHitNoS[0] = std::make_unique<TProfile>(
"QVsNoS0",
";Hit No;Q", 1000, 0., 1000., 0., 2000.);
120 fChargeVsHitNoS[1] = std::make_unique<TProfile>(
"QVsNoS1",
";Hit No;Q", 1000, 0., 1000., 0., 2000.);
121 fChargeVsHitNoS[2] = std::make_unique<TProfile>(
"QVsNoS2",
";Hit No;Q", 1000, 0., 1000., 0., 2000.);
123 fBadWPulseHeight = std::make_unique<TH1D>(
"BWPulseHeight",
"PH (ADC)", 300, 0., 150.);
124 fBadWPulseHVsWidth = std::make_unique<TH2D>(
"BWPHVsWidth",
";PH;Width", 100, 0., 100., 100, 0., 10.);
127 fSPHvsIdx[0] = std::make_unique<TH2D>(
"SPHVsIdx0",
";PH;Idx", 30, 0., 30., 100, 0., 100.);
128 fSPHvsIdx[1] = std::make_unique<TH2D>(
"SPHVsIdx1",
";PH;Idx", 30, 0., 30., 100, 0., 100.);
129 fSPHvsIdx[2] = std::make_unique<TH2D>(
"SPHVsIdx2",
";PH;Idx", 30, 0., 30., 100, 0., 100.);
131 fSWidVsIdx[0] = std::make_unique<TH2D>(
"SWidsIdx0",
";Width;Idx", 30, 0., 30., 100, 0., 10.);
132 fSWidVsIdx[1] = std::make_unique<TH2D>(
"SWidsIdx1",
";Width;Idx", 30, 0., 30., 100, 0., 10.);
133 fSWidVsIdx[2] = std::make_unique<TH2D>(
"SWidsIdx2",
";Width;Idx", 30, 0., 30., 100, 0., 10.);
135 f1PPHvsWid[0] = std::make_unique<TH2D>(
"1PPHVsWid0",
";PH;Width", 100, 0., 100., 100, 0., 20.);
136 f1PPHvsWid[1] = std::make_unique<TH2D>(
"1PPHVsWid1",
";PH;Width", 100, 0., 100., 100, 0., 20.);
137 f1PPHvsWid[2] = std::make_unique<TH2D>(
"1PPHVsWid2",
";PH;Width", 100, 0., 100., 100, 0., 20.);
139 fSPPHvsWid[0] = std::make_unique<TH2D>(
"SPPHVsWid0",
";PH;Width", 100, 0., 100., 100, 0., 20.);
140 fSPPHvsWid[1] = std::make_unique<TH2D>(
"SPPHVsWid1",
";PH;Width", 100, 0., 100., 100, 0., 20.);
141 fSPPHvsWid[2] = std::make_unique<TH2D>(
"SPPHVsWid2",
";PH;Width", 100, 0., 100., 100, 0., 20.);
143 fSOPHvsWid[0] = std::make_unique<TH2D>(
"SOPHVsWid0",
";PH;Width", 100, 0., 100., 100, 0., 20.);
144 fSOPHvsWid[1] = std::make_unique<TH2D>(
"SOPHVsWid1",
";PH;Width", 100, 0., 100., 100, 0., 20.);
145 fSOPHvsWid[2] = std::make_unique<TH2D>(
"SOPHVsWid2",
";PH;Width", 100, 0., 100., 100, 0., 20.);
147 fPHRatVsIdx[0] = std::make_unique<TH2D>(
"PHRatVsIdx0",
";PHRat;Idx", 30, 0., 30., 51, 0., 1.02);
148 fPHRatVsIdx[1] = std::make_unique<TH2D>(
"PHRatVsIdx1",
";PHRat;Idx", 30, 0., 30., 51, 0., 1.02);
149 fPHRatVsIdx[2] = std::make_unique<TH2D>(
"PHRatVsIdx2",
";PHRat;Idx", 30, 0., 30., 51, 0., 1.02);
152 for(
int idx = 0; idx < 3; idx++)
193 size_t longTrackLen(0);
195 for(
const auto& trackHitVecMapItr : trackPlaneHitMap)
199 for(
const auto& planeHitPair : trackHitVecMapItr.second)
203 if (planeHitPair.first == 2 && planeHitPair.second.size() > numHits) numHits = planeHitPair.second.size();
206 if (numHits > longTrackLen)
208 longTrackID = trackHitVecMapItr.first;
209 longTrackLen = numHits;
213 if (longTrackLen > 0)
215 for(
const auto& planeHitPair : trackPlaneHitMap.find(longTrackID)->second)
219 for(
const auto&
hit : planeHitPair.second)
221 if (
hit.Multiplicity() < 2)
fChargeVsHitNoS[planeHitPair.first]->Fill(
float(hitNo)+0.5, std::min(
float(1999.),
hit.Integral()), 1.);
222 fPulseHVsHitNo[planeHitPair.first]->Fill(
float(hitNo)+0.5, std::min(
float(99.9),
hit.PeakAmplitude()), 1.);
223 fChargeVsHitNo[planeHitPair.first]->Fill(
float(hitNo)+0.5, std::min(
float(1999.),
hit.Integral()), 1.);
235 size_t nHitsPerPlane[] = {0,0,0};
238 std::vector<const recob::Hit*> hitSnippetVec;
241 for(
const auto&
hit : hitVec)
245 float chi2DOF = std::min(
hit.GoodnessOfFit(),float(249.8));
246 int numDOF =
hit.DegreesOfFreedom();
247 int hitMult =
hit.Multiplicity();
248 float peakTime =
hit.PeakTime();
249 float charge =
hit.Integral();
250 float sumADC =
hit.SummedADC();
251 float hitPH = std::min(
hit.PeakAmplitude(),float(249.8));
252 float hitSigma =
hit.RMS();
254 size_t plane = wireID.
Plane;
255 size_t wire = wireID.
Wire;
263 std::cout <<
"Hit plane: " << plane <<
", wire: " << wire <<
", T: " << peakTime <<
", PH: " << hitPH <<
", charge: " << charge <<
", sumADC: " << sumADC << std::endl;
266 nHitsPerPlane[plane]++;
275 fFitWidth[plane]->Fill(std::min(
float(19.99),hitSigma), 1.);
284 fPulseHVsWidth[plane]->Fill(std::min(
float(99.9),hitPH), std::min(
float(19.99),hitSigma), 1.);
285 fPulseHVsCharge[plane]->Fill(std::min(
float(99.9),hitPH), std::min(
float(1999.),charge), 1.);
287 if (plane == 2 && hitPH < 5 && hitSigma < 2.2)
289 std::cout <<
"++> plane: " << plane <<
", wire: " << wire <<
", peakTime: " << peakTime <<
", ph: " << hitPH <<
", w: " << hitSigma << std::endl;
292 fBadWPulseHVsWidth->Fill(std::min(
float(99.9),hitPH), std::min(
float(19.99),hitSigma), 1.);
300 if (!hitSnippetVec.empty() && hitSnippetVec.back()->LocalIndex() >=
hit.LocalIndex())
303 if (hitSnippetVec.size() > 1)
306 std::sort(hitSnippetVec.begin(),hitSnippetVec.end(),[](
const auto*
left,
const auto*
right){
return left->PeakAmplitude() >
right->PeakAmplitude();});
308 float maxPulseHeight = hitSnippetVec.front()->PeakAmplitude();
310 for(
size_t idx = 0; idx < hitSnippetVec.size(); idx++)
312 float pulseHeight = hitSnippetVec.at(idx)->PeakAmplitude();
313 float pulseWid = hitSnippetVec.at(idx)->RMS();
314 float pulseHeightRatio = pulseHeight / maxPulseHeight;
316 size_t plane = hitSnippetVec.at(idx)->WireID().Plane;
318 fSPHvsIdx[plane]->Fill(idx, std::min(
float(99.9),pulseHeight), 1.);
319 fSWidVsIdx[plane]->Fill(idx, std::min(
float(19.99),pulseWid), 1.);
320 fPHRatVsIdx[plane]->Fill(idx, pulseHeightRatio, 1.);
322 if (idx == 0)
fSPPHvsWid[plane]->Fill(std::min(
float(99.9),pulseHeight), std::min(
float(19.99),pulseWid), 1.);
323 else fSOPHvsWid[plane]->Fill(std::min(
float(99.9),pulseHeight), std::min(
float(19.99),pulseWid), 1.);
328 float pulseHeight = hitSnippetVec.front()->PeakAmplitude();
329 float pulseWid = hitSnippetVec.front()->RMS();
330 size_t plane = hitSnippetVec.front()->WireID().Plane;
332 f1PPHvsWid[plane]->Fill(std::min(
float(99.9),pulseHeight), std::min(
float(19.99),pulseWid), 1.);
335 hitSnippetVec.clear();
338 hitSnippetVec.push_back(&
hit);
348 double normFactor(1./numEvents);
356 for(
int idx = 0; idx < 3; idx++)
std::unique_ptr< TH1D > fHitSumADC[3]
std::unique_ptr< TH1D > fPulseHeightSingle[3]
const geo::GeometryCore * geometry
~HitAnalysisAlg()
Destructor.
std::unique_ptr< TH1D > fChi2DOF[3]
std::unique_ptr< TH1D > fChi2DOFSingle[3]
void endJob(int numEvents)
std::vector< recob::Hit > HitVec
std::unique_ptr< TProfile > fChargeVsHitNoS[3]
std::unique_ptr< TProfile > fPulseHVsHitNo[3]
std::unique_ptr< TH1D > fHitCharge[3]
WireID_t Wire
Index of the wire within its plane.
Scale(size_t pos, T factor) -> Scale< T >
const geo::GeometryCore * fGeometry
pointer to Geometry service
std::unique_ptr< TH1D > fHitsByWire[3]
std::unique_ptr< TH2D > fSWidVsIdx[3]
unsigned int Nwires(unsigned int p, unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wires in the specified plane.
std::unique_ptr< TH1D > fNumDegFree[3]
std::unique_ptr< TH2D > fPulseHVsCharge[3]
std::unique_ptr< TH2D > f1PPHvsWid[3]
std::unique_ptr< TProfile > fChargeVsHitNo[3]
std::unique_ptr< TH1D > fHitMult[3]
std::string fLocalDirName
Fraction for truncated mean.
std::unique_ptr< TH1D > fPulseHeightMulti[3]
std::unique_ptr< TH1D > fBadWPulseHeight
PlaneID_t Plane
Index of the plane within its TPC.
std::unique_ptr< TH1D > fPulseHeight[3]
Description of geometry of one entire detector.
TDirectory * fRootDirectory
std::unique_ptr< TH2D > fSPPHvsWid[3]
std::unique_ptr< TH2D > fSPHvsIdx[3]
std::unique_ptr< TH2D > fPHRatVsIdx[3]
std::unique_ptr< TH1D > fFitWidth[3]
void fillHistograms(const TrackPlaneHitMap &) const
std::unique_ptr< TH1D > fDriftTimes[3]
std::map< int, PlaneHitMap > TrackPlaneHitMap
std::unique_ptr< TH1D > fHitsByTime[3]
std::unique_ptr< TH2D > fNDFVsChi2[3]
HitAnalysisAlg(fhicl::ParameterSet const &pset)
void setup(const geo::GeometryCore &, TDirectory *)
Begin job method.
std::unique_ptr< TH2D > fBadWPulseHVsWidth
std::unique_ptr< TH2D > fSOPHvsWid[3]
std::unique_ptr< TH1D > fBadWHitsByWire
BEGIN_PROLOG could also be cout
std::unique_ptr< TH2D > fPulseHVsWidth[3]
void reconfigure(fhicl::ParameterSet const &pset)