SNSlicerAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SNSlicerAna
3 // Module Type: analyzer
4 // File: SNSlicerAna_module.cc
5 //
6 // Generated at Tue May 23 12:56:03 2017 by Justin Vasel using artmod
7 // from cetpkgsupport v1_10_02.
8 ////////////////////////////////////////////////////////////////////////
9 
10 // ROOT includes
11 #include <TTree.h>
12 #include <TVector3.h>
13 
14 // Framework includes
21 #include "fhiclcpp/ParameterSet.h"
23 
24 // NOvASoft includes
25 #include "Geometry/Geometry.h"
28 #include "MCCheater/BackTracker.h"
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/Cluster.h"
31 #include "Simulation/Particle.h"
32 
33 
34 namespace sn {
35  class SNSlicerAna;
36 }
37 
39 public:
40  explicit SNSlicerAna(fhicl::ParameterSet const & p);
41 
42  SNSlicerAna(SNSlicerAna const &) = delete;
43  SNSlicerAna(SNSlicerAna &&) = delete;
44  SNSlicerAna & operator = (SNSlicerAna const &) = delete;
45  SNSlicerAna & operator = (SNSlicerAna &&) = delete;
46 
47  // Required functions.
48  void analyze(art::Event const & e) override;
49 
50  // Selected optional functions.
51  void beginJob() override;
52 
53 private:
56 
57 
60 
61  TTree* fClusterTree;
62  TTree* fEventTree;
63  TTree* fHitTree;
64  TTree* fHitPairTree;
65 
66  // ART stuff
67  int fEventNumber; /// ART event number
68  int fClusterId; /// Unique ID for each cluster in the event
69  double fClusterTNS; /// Min time of the cluster
70 
71  // Hit metrics
72  int fNHitsMC; /// Number of MC hits in cluster
73  int fNHitsBG; /// Number of non-MC hits in cluster
74  int fAdcMc; /// ADC for MC hits
75  int fAdcBg; /// ADC for BG hits
76  int fNFEB; /// Number of distinct FEBs in the cluster
77 
78  // Cluster metrics
79  double fSumADC; /// Total cluster ADC
80  int fSpanADC; /// Difference between max and min ADCs
81  double fMeanADC; /// Mean cluster ADC
82  double fStdADC; /// Std. Dev. of cluster ADCs
83  double fSumMeV; /// Cluster total energy (MeV)
84  double fSpanTNS; /// Cluster time span (ns)
85  int fSpanPlane; /// Cluster plane span (planes)
86  int fSpanCellX; /// Cluster cell span in X-view (cells)
87  int fSpanCellY; /// Cluster cell span in Y-view (cells)
88  int fSpanViews; /// Cluster view span (1 or 2)
89  float fPairwiseADC; /// Average of hit-pairwise ADC differences in a cluster
90  float fPairwiseTNS; /// Average of hit-pairwise TNS differences in a cluster
91 
92  // Means
93  double fMeanTNS; /// Cluster time mean (ns)
94  double fMeanX; /// Cluster X mean
95  double fMeanY; /// Cluster Y mean
96  double fMeanZ; /// Cluster Z mean
97 
98  float fClusterEfficiency; /// Number of MC hits (in all clusters / total in event)
99 
100  float fPurity; /// Cluster Purity (Nsig / Ntotal)
101  int fAdcDiff; /// Difference in ADC between pair-wise hits
104  float fXdiff;
105 };
106 
107 
109  EDAnalyzer(p),
110  fClusterLabel(p.get<std::string>("ClusterLabel")),
111  fCellHitLabel(p.get<std::string>("CellHitLabel")),
112  fEventNumber(0),
113  fClusterId(0),
114  fClusterTNS(0),
115  fNHitsMC(0),
116  fNHitsBG(0),
117  fAdcMc(0),
118  fAdcBg(0),
119  fNFEB(0),
120  fSumADC(0),
121  fSpanADC(0),
122  fMeanADC(0),
123  fStdADC(0),
124  fSumMeV(0),
125  fSpanTNS(0),
126  fSpanPlane(0),
127  fSpanCellX(0),
128  fSpanCellY(0),
129  fSpanViews(0),
130  fPairwiseADC(0),
131  fPairwiseTNS(0),
132  fMeanTNS(0),
133  fMeanX(0),
134  fMeanY(0),
135  fMeanZ(0),
137  fPurity(0),
138  fAdcDiff(0),
139  fPlaneDiff(0),
140  fCellDiff(0),
141  fXdiff(0)
142 {
143 }
144 
145 
147 {
148  fEventNumber = e.id().event();
149 
150  // Get the clusters from the ART event
152  e.getByLabel(fClusterLabel, clusters);
153 
154  int fClusterHitsMC = 0;
155 
156  for (size_t iCluster = 0; iCluster < clusters->size(); ++iCluster) {
157  rb::Cluster aCluster = clusters->at(iCluster);
158  fClusterId = iCluster;
159 
160  if (aCluster.NCell()==0) continue;
161 
162  fClusterTNS = aCluster.MinTNS();
163 
164  fSumADC = aCluster.TotalADC();
165  fSumMeV = aCluster.TotalGeV() * 1.0e3;
166  fMeanADC = aCluster.TotalADC() / (float)aCluster.NCell();
167  fSpanTNS = aCluster.ExtentTNS();
168  fSpanPlane = aCluster.ExtentPlane();
169  fSpanCellX = (aCluster.NXCell()) ? aCluster.ExtentCell(geo::kX) : 0;
170  fSpanCellY = (aCluster.NYCell()) ? aCluster.ExtentCell(geo::kY) : 0;
171 
172  fSpanViews = (fSpanCellX>0 && fSpanCellY>0) ? 2 : 1;
173 
174  fMeanTNS = aCluster.MeanTNS();
175  fMeanX = aCluster.MeanX();
176  fMeanY = aCluster.MeanY();
177  fMeanZ = aCluster.MeanZ();
178 
179 
180  // Loop over the contained cell hits
181  std::map<int, int> mapInteractionCounts;
182  fNHitsMC = 0;
183  fNHitsBG = 0;
184  std::vector<unsigned int> fFEBs;
185 
186  int fMinADC=0;
187  int fMaxADC=0;
188  float fVarADC=0;
189  art::PtrVector<rb::CellHit> allHits = aCluster.AllCells();
190  for (art::Ptr<rb::CellHit> aHit : allHits) {
191  fAdcMc = 0;
192  fAdcBg = 0;
193  (aHit->IsMC()) ? fNHitsMC++ : fNHitsBG++;
194  (aHit->IsMC()) ? fAdcMc = aHit->ADC() : fAdcBg = aHit->ADC();
195 
196  if (fMinADC==0 || aHit->ADC() < fMinADC) fMinADC = aHit->ADC();
197  if (fMaxADC==0 || aHit->ADC() > fMaxADC) fMaxADC = aHit->ADC();
198 
199  fVarADC += pow(aHit->ADC() - fMeanADC, 2) / (float)aCluster.NCell();
200 
201  unsigned int fFEB = aHit->DaqChannel() & 0xFFFFFF00;
202  if (std::find(fFEBs.begin(), fFEBs.end(), fFEB) == fFEBs.end()) fFEBs.push_back(fFEB);
203 
204  // fHitTree->Fill();
205  }
206 
207  fStdADC = sqrt(fVarADC);
208 
209  int n = 0;
210  float totalPairwiseADC = 0;
211  float totalPairwiseTNS = 0;
212  for (size_t i=0; i<allHits.size()-1; ++i) {
213  for (size_t j=i+1; j<allHits.size(); ++j) {
214  ++n;
215 
216  totalPairwiseADC = std::abs(allHits[i]->ADC() - allHits[j]->ADC());
217  totalPairwiseTNS = std::abs(allHits[i]->TNS() - allHits[j]->TNS());
218  fAdcDiff = std::abs(allHits[i]->ADC() - allHits[j]->ADC());
219  fPlaneDiff = std::abs(allHits[i]->Plane() - allHits[j]->Plane());
220  fCellDiff = std::abs(allHits[i]->Cell() - allHits[j]->Cell());
222 
223  TVector3 x1, x2;
224  fGeom->Plane(allHits[i]->Plane())->Cell(allHits[i]->Cell())->GetCenter(x1);
225  fGeom->Plane(allHits[j]->Plane())->Cell(allHits[j]->Cell())->GetCenter(x2);
226 
227  fXdiff = (x1 - x2).Mag();
228 
229  fHitPairTree->Fill();
230  }
231  }
232  fPairwiseADC = totalPairwiseADC / (float)n;
233  fPairwiseTNS = totalPairwiseTNS / (float)n;
234 
235  fSpanADC = fMaxADC - fMinADC;
236 
237  fNFEB = fFEBs.size();
238  fClusterHitsMC += fNHitsMC;
239 
240  fClusterTree->Fill();
241  }
242 
243  // Get the hits from the ART event
245  e.getByLabel(fCellHitLabel, hits);
246 
247  int fEventHitsMC = 0;
248  for (rb::CellHit aHit : *hits) {
249  if (aHit.IsMC()) ++fEventHitsMC;
250  }
251 
252  fClusterEfficiency = (float)fClusterHitsMC / (float)fEventHitsMC;
253 
254  fEventTree->Fill();
255 }
256 
257 
259 {
260  fClusterTree = tfs->make<TTree>("ClusterTree", "ClusterTree");
261  fClusterTree->Branch("EventNumber", &fEventNumber);
262  fClusterTree->Branch("ClusterID", &fClusterId);
263  fClusterTree->Branch("ClusterTNS", &fClusterTNS);
264  fClusterTree->Branch("NHitsMC", &fNHitsMC);
265  fClusterTree->Branch("NHitsBG", &fNHitsBG);
266  fClusterTree->Branch("NFEB", &fNFEB);
267  fClusterTree->Branch("SumADC", &fSumADC);
268  fClusterTree->Branch("SpanADC", &fSpanADC);
269  fClusterTree->Branch("MeanADC", &fMeanADC);
270  fClusterTree->Branch("StdADC", &fStdADC);
271  fClusterTree->Branch("SumMeV", &fSumMeV);
272  fClusterTree->Branch("PairwiseADC", &fPairwiseADC);
273  fClusterTree->Branch("PairwiseTNS", &fPairwiseTNS);
274  fClusterTree->Branch("SpanTNS", &fSpanTNS);
275  fClusterTree->Branch("SpanPlane", &fSpanPlane);
276  fClusterTree->Branch("SpanCellX", &fSpanCellX);
277  fClusterTree->Branch("SpanCellY", &fSpanCellY);
278  fClusterTree->Branch("SpanViews", &fSpanViews);
279  fClusterTree->Branch("MeanTNS", &fMeanTNS);
280  fClusterTree->Branch("MeanX", &fMeanX);
281  fClusterTree->Branch("MeanY", &fMeanY);
282  fClusterTree->Branch("MeanZ", &fMeanZ);
283 
284  fEventTree = tfs->make<TTree>("EventTree", "EventTree");
285  fEventTree->Branch("EventNumber", &fEventNumber);
286  fEventTree->Branch("ClusterEfficiency", &fClusterEfficiency);
287 
288  fHitTree = tfs->make<TTree>("HitTree", "HitTree");
289  fHitTree->Branch("EventNumber", &fEventNumber);
290  fHitTree->Branch("AdcMc", &fAdcMc);
291  fHitTree->Branch("AdcBg" , &fAdcBg);
292 
293  fHitPairTree = tfs->make<TTree>("HitPairTree", "HitPairTree");
294  fHitPairTree->Branch("EventNumber", &fEventNumber);
295  fHitPairTree->Branch("ClusterID", &fClusterId);
296  fHitPairTree->Branch("ClusterTNS", &fClusterTNS);
297  fHitPairTree->Branch("Purity", &fPurity);
298  fHitPairTree->Branch("PlaneDiff", &fPlaneDiff);
299  fHitPairTree->Branch("CellDiff", &fCellDiff);
300  fHitPairTree->Branch("ADCdiff", &fAdcDiff);
301  fHitPairTree->Branch("Xdiff", &fXdiff);
302 }
303 
back track the reconstruction to the simulation
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
int fAdcDiff
Cluster Purity (Nsig / Ntotal)
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
void analyze(art::Event const &e) override
double fMeanTNS
Average of hit-pairwise TNS differences in a cluster.
float fPairwiseADC
Cluster view span (1 or 2)
double fSumMeV
Std. Dev. of cluster ADCs.
art::ServiceHandle< art::TFileService > tfs
Float_t x1[n_points_granero]
Definition: compare.C:5
int fPlaneDiff
Difference in ADC between pair-wise hits.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:233
int fNFEB
ADC for BG hits.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
SNSlicerAna(fhicl::ParameterSet const &p)
int fSpanViews
Cluster cell span in Y-view (cells)
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
double fMeanZ
Cluster Y mean.
unsigned int ExtentCell(geo::View_t view) const
Definition: Cluster.cxx:570
double MeanX(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:231
double MinTNS() const
Definition: Cluster.cxx:482
double TotalADC() const
Sum of the ADC of all the contained hits.
Definition: Cluster.cxx:360
void hits()
Definition: readHits.C:15
double MeanY(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:232
int fSpanCellY
Cluster cell span in X-view (cells)
double fMeanY
Cluster X mean.
int fNHitsMC
Min time of the cluster.
double fMeanX
Cluster time mean (ns)
SNSlicerAna & operator=(SNSlicerAna const &)=delete
Remove hits from hot and cold channels.
const double j
Definition: BetheBloch.cxx:29
int fSpanCellX
Cluster plane span (planes)
int fClusterId
ART event number.
float fPairwiseTNS
Average of hit-pairwise ADC differences in a cluster.
std::string fCellHitLabel
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double fSumADC
Number of distinct FEBs in the cluster.
float fPurity
Number of MC hits (in all clusters / total in event)
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
int fNHitsBG
Number of MC hits in cluster.
EventNumber_t event() const
Definition: EventID.h:116
int fAdcMc
Number of non-MC hits in cluster.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
art::ServiceHandle< geo::Geometry > fGeom
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double fClusterTNS
Unique ID for each cluster in the event.
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
float Mag() const
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
int fSpanPlane
Cluster time span (ns)
double fMeanADC
Difference between max and min ADCs.
double fSpanTNS
Cluster total energy (MeV)
Float_t e
Definition: plot.C:35
int fAdcBg
ADC for MC hits.
double ExtentTNS() const
Definition: Cluster.h:252
int fSpanADC
Total cluster ADC.
std::string fClusterLabel
double fStdADC
Mean cluster ADC.
float fClusterEfficiency
Cluster Z mean.
EventID id() const
Definition: Event.h:56
void beginJob() override
Encapsulate the geometry of one entire detector (near, far, ndos)