ClusterAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: ClusterAna
3 // Module Type: analyzer
4 // File: ClusterAna_module.cc
5 //
6 // Generated at Tue Apr 2 09:40:02 2013 by Jan Zirnstein using artmod
7 // from art v1_02_06.
8 ////////////////////////////////////////////////////////////////////////
9 
10 // Framework includes
16 
17 // NOvADDT includes
23 
24 // ROOT includes
25 #include "TH1F.h"
26 
27 #include <algorithm>
28 #include <cmath>
29 #include <vector>
30 
31 namespace novaddt {
32  class ClusterAna;
33 }
34 
36 public:
37  explicit ClusterAna(fhicl::ParameterSet const & p);
38  virtual ~ClusterAna();
39 
40  void analyze(art::Event const & e) override;
41 
42  void beginJob() override;
43  void endJob() override;
44 
45 private:
46 
51  TH1F* fWidthXHist;
52  TH1F* fWidthYHist;
57 
59 
60 };
61 
62 
64  EDAnalyzer(p),
65  fClusterModule(p.get<std::string>("ClusterModule"))
66 {
67 }
68 
70 {
71  // Clean up dynamic memory and other resources here.
72 }
73 
75 {
77  e.getByLabel(fClusterModule, clusters);
78 
79  // put in a mf::LogWarning why this statement returns if true
80  if(clusters->empty()){
81  mf::LogWarning ("No Clusters")<<"No Clusters in the input file";
82  return;
83  }
84 
85  // Fill Number of Clusters in Event Histogram
86  fNClustersHist->Fill(clusters->size());
87 
88  // Grab the Hitlists associated with the clusters
90 
91  //std::cout<< "Sanity check, there are: "<<clusters->size()<<" clusters and " <<hls.size()<<" hitlists."<<std::endl;
92 
93  unsigned int i(0);
94  for(auto const& clust: *clusters){
95  fExtentPlanesHist->Fill(clust.Extent());
96  fWidthXHist->Fill(clust.WidthX());
97  fWidthYHist->Fill(clust.WidthY());
98  fClusterEnergyHist->Fill(clust.TotalADC());
99  fClusterEnDenHist->Fill(clust.TotalADC()/(1.0*clust.Volume()));
100 
101  art::Ptr<novaddt::HitList> hl_ptr(hls.at(i));
102  ++i;
103 
104  novaddt::HitList hlist(*hl_ptr);
105 
106  fNHitsPerClusterHist->Fill(hl_ptr->size());
107  //std::cout<<"Hits in cluster: "<<hlist.size()<<std::endl;
108 
109  // Split the hits in x and y to do 2D fitting
110  auto bound = std::partition(hlist.begin(),hlist.end(),CompareDAQHit<View>());
111  //std::cout<<"Boundary is at: "<<bound-hlist.begin()<<std::endl;
112 
113  unsigned int xsize = bound-hlist.begin();
114  unsigned int ysize = hlist.end()-bound;
115 
116  //std::cout<<"X Size: "<<xsize<<", Y size: "<<ysize<<std::endl;
117 
118  // Smarter way of vector managment
119 
120  std::vector<double> XXpoints;
121  std::vector<double> XZpoints;
122  std::vector<double> YYpoints;
123  std::vector<double> YZpoints;
124 
125  XXpoints.reserve(xsize);
126  XZpoints.reserve(xsize);
127  YYpoints.reserve(ysize);
128  YZpoints.reserve(ysize);
129 
130  {
131  auto it = hlist.begin();
132  for(it=hlist.begin();it!=bound;++it){
133  XXpoints.emplace_back((*it).Cell().val);
134  XZpoints.emplace_back((*it).Plane().val);
135  }
136 
137  for(it=bound;it!=hlist.end();++it){
138  YYpoints.emplace_back((*it).Cell().val);
139  YZpoints.emplace_back((*it).Plane().val);
140  }
141  }
142 
143  // Fit to a straight line
144  double xslope;
145  double xint;
146  unsigned int xdof = xsize-1;
147  double xsigm=0.0;
148  double pred;
149 
150  if(xsize > 1){
151  util::LinFitUnweighted(XXpoints, XZpoints, xslope, xint);
152  //std::cout<<"XView fit done: y = "<<xslope<<"*x + "<<xint<<std::endl;
153 
154  // Evaluate goodness of fit
155  for(unsigned int xi=0;xi<xsize;++xi){
156  pred = xslope*XXpoints[xi] + xint;
157  xsigm += sqrt(util::sqr(XZpoints[xi]-pred));
158  }
159 
160  fXViewTrackGood->Fill(xsigm/xdof);
161  }
162 
163  else{
164  fXViewTrackGood->Fill(0.0);
165  xdof = 1;
166  }
167  //std::cout<<"Goodness of XView fit: "<<xsigm/xdof<<std::endl;
168 
169  double yslope;
170  double yint;
171  unsigned int ydof = ysize-1;
172  double ysigm=0.0;
173 
174  if(ysize > 1){
175  util::LinFitUnweighted(YYpoints, YZpoints, yslope, yint);
176  //std::cout<<"YView fit done: y = "<<yslope<<"*x + "<<yint<<std::endl;
177 
178  for(unsigned int yi=0;yi<ysize;++yi){
179  pred = yslope*YYpoints[yi] + yint;
180  ysigm += sqrt(util::sqr(YZpoints[yi]-pred));
181  }
182 
183  fYViewTrackGood->Fill(ysigm/ydof);
184  }
185 
186  else{
187  fYViewTrackGood->Fill(0.0);
188  ydof = 1;
189  }
190  //std::cout<<"Goodness of YView fit: "<<ysigm/ydof<<std::endl;
191 
192  fCombViewTrackGood->Fill(sqrt(util::sqr(xsigm/xdof)+util::sqr(ysigm/ydof)));
193 
194  //std::cout<<"Analysis of Cluster finished\n\n";
195 
196  } // end of cluster
197 
198 } // end of analyze function
199 
201 {
203 
204  fNClustersHist = f->make<TH1F>("fNClusters",";NClusters;Events",100,0.,100.);
205  fNHitsPerClusterHist = f->make<TH1F>("fNHitsPerCluster",";NHits;Clusters",100,0.,100.);
206  fExtentPlanesHist = f->make<TH1F>("fExtentPlanes",";NPlanes;Clusters",896,0.,896.);
207  fWidthXHist = f->make<TH1F>("fWidthX",";NCells;Clusters",384,0.,384.);
208  fWidthYHist = f->make<TH1F>("fWidthY",";NCells;Clusters",384,0.,384.);
209  fClusterEnergyHist = f->make<TH1F>("fClusterEnergy",";TotADC;Clusters",1000,0.,100000.);
210  fClusterEnDenHist = f->make<TH1F>("fClusterEnergyDensity",";ADC per plane*cell^2;Clusters",10000,0.,0.5);
211  fXViewTrackGood = f->make<TH1F>("fXViewTrackGood",";Goodness of Lin Fit;Tracks",1000,0.,100.);
212  fYViewTrackGood = f->make<TH1F>("fYViewTrackGood",";Goodness of Lin Fit;Tracks",1000,0.,100.);
213  fCombViewTrackGood = f->make<TH1F>("fCombViewTrackGood",";Goodness of Lin Fit;Tracks",1000,0.,100.);
214 
215 }
216 
218 {
219  // Implementation of optional member function here.
220 }
221 
ClusterAna(fhicl::ParameterSet const &p)
set< int >::iterator it
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
DEFINE_ART_MODULE(TestTMapFile)
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
void endJob() override
void beginJob() override
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void LinFitUnweighted(const std::vector< double > &x, const std::vector< double > &y, double &m, double &c)
Simplified version of LinFit.
Definition: MathUtil.cxx:8
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Float_t e
Definition: plot.C:35
Definition: fwd.h:28
void analyze(art::Event const &e) override