HitEfficiency_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: HitEfficiency
3 // Module Type: Analyzer
4 // File: HitEfficiency_module.cc
5 //
6 ////////////////////////////////////////////////////////////////////////
7 
8 // framework includes....
15 #include "fhiclcpp/ParameterSet.h"
19 
20 // novasoft stuff
21 #include "RecoBase/Track.h"
22 #include "RecoBase/CellHit.h"
23 #include "Calibrator/Calibrator.h"
24 #include "Geometry/Geometry.h"
25 #include "Geometry/LiveGeometry.h"
27 #include "Utilities/func/MathUtil.h"
29 #include "RunHistory/RunHistory.h"
31 
32 
33 // C++
34 #include <memory>
35 #include <string>
36 #include <iostream>
37 #include <vector>
38 
39 // ROOT
40 #include "TVector3.h"
41 #include "TH2F.h"
42 
43 
44 namespace calib {
45 
46  class HitEfficiency : public art::EDAnalyzer {
47  public:
48  explicit HitEfficiency(fhicl::ParameterSet const & p);
49  virtual ~HitEfficiency();
50 
51  void analyze(art::Event const& e);
52  void reconfigure(fhicl::ParameterSet const & p);
53  void beginJob();
54  void beginRun(art::Run const& r);
55  void endJob();
56 
57  private:
58 
60  double fNTrkHits;
61  double fNPlanes;
62  double fMaxPlane;
63  double fDEdge;
64 
65  unsigned int fplanemax; // maximum plane to calculate efficiencies on
66  double zmaxinstr; // maximum instrumented z position
67  double zmininstr; // minimum instrumented z position
68 
69  //int fDet; // Detector ID
71 
72  TH2F* hitsVwVlenx;
73  TH2F* cellsVwVlenx;
74 
75  TH2F* hitsVwVleny;
76  TH2F* cellsVwVleny;
77 
80 
83 
84  // efficiency plots
85  TH2F* hitsVwcell;
86  TH2F* cellsVwcell;
87 
88  // plots after applying track quality cuts plotted in relation to distance from apd, not w
91 
92  // plots after applying path length cut
95 
98 
99 
100  // Histograms to keep track of the number of events used
101  TH1F* NEvents;
102  int fnevents;
103 
104  // define length cuts
105  double minpath = 6.0;
106  double maxpath = 7.0;
107 
108  };
109 
110 
112  : EDAnalyzer(p)
113  {
114  this->reconfigure(p);
115  }
116 
118  {
119  fTrackLabel = pset.get< std::string >("TrackLabel");
120  fNTrkHits = pset.get< double >("NTrkHits");
121  fNPlanes = pset.get< double >("NPlanes");
122  fMaxPlane = pset.get< double >("MaxPlane");
123  fDEdge = pset.get< double >("DEdge");
124  }
125 
126 
128  {
129  // Clean up dynamic memory and other resources here.
130  }
131 
133  {
135 
136  fplanemax = 896;
137 
138  fnevents = 0;
139 
140  zmaxinstr = 0.0; // maximum instrumented z position
141  zmininstr = 0.0;
142 
143  // Define my histograms binning here
144  double nwbins = 32;
145  double wmin = -800;
146  double wmax = 800;
147 
148  double napdbins = nwbins;
149  double apdmin = 0.0;
150  double apdmax = 1600;
151 
152  double nlbins = 25;
153  double lmin = 0.0;
154  double lmax = 10.0;
155 
156  int ncellbins = 384*fplanemax+1;
157  double cellmin = -0.5;
158  double cellmax = 384*(double)fplanemax+0.5;
159 
160  // efficiency plots
161  hitsVwVlenx = tfs->make<TH2F>("hitsVwVlenx",";W (cm);Length (cm)",
162  nwbins,wmin,wmax,nlbins,lmin,lmax);
163  cellsVwVlenx = tfs->make<TH2F>("cellsVwVlenx",";W (cm);Length (cm)",
164  nwbins,wmin,wmax,nlbins,lmin,lmax);
165 
166 
167  hitsVwVleny = tfs->make<TH2F>("hitsVwVleny",";W (cm);Length (cm)",
168  nwbins,wmin,wmax,nlbins,lmin,lmax);
169  cellsVwVleny = tfs->make<TH2F>("cellsVwVleny",";W (cm);Length (cm)",
170  nwbins,wmin,wmax,nlbins,lmin,lmax);
171 
172 
173  hitsVapddistVlenx = tfs->make<TH2F>("hitsVapddistVlenx",";Distance from APD (cm);Length (cm)",
174  napdbins,apdmin,apdmax,nlbins,lmin,lmax);
175  cellsVapddistVlenx = tfs->make<TH2F>("cellsVapddistVlenx",";Distance from APD (cm);Length (cm)",
176  napdbins,apdmin,apdmax,nlbins,lmin,lmax);
177 
178 
179  hitsVapddistVleny = tfs->make<TH2F>("hitsVapddisVleny",";Distance from APD (cm);Length (cm)",
180  napdbins,apdmin,apdmax,nlbins,lmin,lmax);
181  cellsVapddistVleny = tfs->make<TH2F>("cellsVapddisVleny",";Distance from APD (cm);Length (cm)",
182  napdbins,apdmin,apdmax,nlbins,lmin,lmax);
183 
184 
185  // efficiency plots as a function of w
186  hitsVwcell = tfs->make<TH2F>("hitsVwcell",";Cell;W (cm)",
187  ncellbins,cellmin,cellmax,nwbins,wmin,wmax);
188  cellsVwcell = tfs->make<TH2F>("cellsVwcell",";Cell;W (cm)",
189  ncellbins,cellmin,cellmax,nwbins,wmin,wmax);
190 
191 
192  // plots as a function of apd distance
193  hitsVapddistcell = tfs->make<TH2F>("hitsVapddistcell",";Cell;Distance from APD (cm)",
194  ncellbins,cellmin,cellmax,napdbins,apdmin,apdmax);
195  cellsVapddistcell = tfs->make<TH2F>("cellsVapddistcell",";Cell;Distance from APD (cm)",
196  ncellbins,cellmin,cellmax,napdbins,apdmin,apdmax);
197 
198 
199  // plots after applying path length cut
200  hitsVwcelllencut = tfs->make<TH2F>("hitsVwcelllencut",";Cell;W (cm)",
201  ncellbins,cellmin,cellmax,nwbins,wmin,wmax);
202  cellsVwcelllencut = tfs->make<TH2F>("cellsVwcelllencut",";Cell;W (cm)",
203  ncellbins,cellmin,cellmax,nwbins,wmin,wmax);
204 
205 
206  hitsVapddistcelllencut = tfs->make<TH2F>("hitsVapddistcelllencut",";Cell;Distance from APD (cm)",
207  ncellbins,cellmin,cellmax,napdbins,apdmin,apdmax);
208  cellsVapddistcelllencut = tfs->make<TH2F>("cellsVapddistcelllencut",";Cell;Distance from APD (cm)",
209  ncellbins,cellmin,cellmax,napdbins,apdmin,apdmax);
210 
211  NEvents = tfs->make<TH1F>("NEvents",";NEvents;",1,-0.5,0.5);
212 
213 
214  }
215 
217  {
218  // make a geometry service handle
220  int planemin = 0;
221  int planemax = geom->NPlanes()-1;
222  if(fMaxPlane >= planemin && fMaxPlane <= planemax){
224  }
225  else{ fplanemax = planemax; }
226 
228  double zfront1,zback1,zfront2,zback2;
229  int ndets = lgeom->InstrumentedDetEnds(zfront1,zback1,zfront2,zback2);
230 
231  // take the most downstream zfront and zback
232  if(ndets > 0){
233  zmaxinstr = zback1;
234  zmininstr = zfront1;
235  }
236 
237  // leaving this disabled for now
238 // fDet = geom->DetId();
239 // // get the run history
240 // run = r.run();
241 // nova::dbi::RunHistory runhist(run,fDet);
242 // frunhist = runhist;
243 // frunhist.LoadDAQRunHistory();
244 // frunhist.LoadPixelInfo();
245 
246  }
247 
249  {
250 
251  ++fnevents;
252 
253  // make a geometry service handle
255  // make a cmap service handle
257 
258 
259  //get the cosmic tracks associated with each event
261  e.getByLabel(fTrackLabel,trackcol);
263  for(unsigned int i = 0; i<trackcol->size(); ++i){
264  art::Ptr<rb::Track>track(trackcol,i);
265  tracks.push_back(track);
266  }
267 
268  // loop over the tracks
269  for(unsigned int itr = 0; itr < tracks.size(); ++itr){
270  art::Ptr<rb::Track> track = tracks[itr];
271 
272  if(track->View() != geo::kXorY){ continue; }
273 
274  // get basic track information here to store in tree
275  int nxhit = track->NXCell();
276  int nyhit = track->NYCell();
277  if(nxhit < fNTrkHits || nyhit < fNTrkHits){ continue; }
278 
279  TVector3 dir = track->Dir();
280  double dx = dir.X();
281  double dy = dir.Y();
282  double dz = dir.Z();
283  if( dx > 1.0 || dx < -1.0 ||
284  dy > 1.0 || dy < -1.0 ||
285  dz > 1.0 || dz < -1.0){ continue; }
286 
287  // find the minimum, maximum planes
288  unsigned int minplane = track->MinPlane();
289  unsigned int maxplane = track->MaxPlane();
290 
291  if((maxplane - minplane) < fNPlanes){ continue; }
292  if(maxplane > fplanemax || minplane > fplanemax){ continue; }
293 
294  // get the track start/stop point info
295  TVector3 start = track->Start();
296  TVector3 stop = track->Stop();
297  double ix = start.X();
298  double iy = start.Y();
299  double iz = start.Z();
300  double fx = stop.X();
301  double fy = stop.Y();
302  double fz = stop.Z();
303 
304  double zmin = zmininstr+fDEdge;
305  double zmax = zmaxinstr-fDEdge;
306  double xmax = geom->DetHalfWidth() - fDEdge;
307  double ymax = geom->DetHalfHeight() - fDEdge;
308 
309  // through going track cut
310  if(fabs(ix) < xmax && fabs(iy) < ymax && iz > zmin && iz < zmax){ continue; } //start
311  if(fabs(fx) < xmax && fabs(fy) < ymax && fz > zmin && fz < zmax){ continue; } // stop
312 
313  // get the hits that belong to this track
314  art::PtrVector<rb::CellHit> trkhits = track->AllCells();
315 
316  // sort the hits by plane
317  rb::SortByPlane(trkhits);
318 
319  // get the list of cells that the track passes through
320  std::vector<geo::OfflineChan> allcellsontrack, ycellsontrack;
321  geom->CountCellsOnLine(start,stop,allcellsontrack,ycellsontrack);
322 
323  for(unsigned int i = 0; i < ycellsontrack.size(); ++i){
324  allcellsontrack.push_back(ycellsontrack[i]);
325  }
326 
327  // make a list of tricells
328  std::vector<geo::OfflineChan> tricells;
329  for(unsigned int icell = 0; icell < allcellsontrack.size(); ++icell){
330  // check if the z position of the cell falls in the instrumented region
331  double zcent[3];
332  geom->Plane(allcellsontrack[icell].Plane())->Cell(allcellsontrack[icell].Cell())->GetCenter(zcent);
333  if(zcent[2] < zmininstr || zcent[2] > zmaxinstr){ continue; }
334  // find the maximum and minimum cell hit on this plane
335  int plane = allcellsontrack[icell].Plane();
336  int cellmin = allcellsontrack[icell].Cell()-1;
337  int cellmax = allcellsontrack[icell].Cell()+1;
338  bool topok = false;
339  bool botok = false;
340  for(unsigned int ihit = 0; ihit < trkhits.size(); ++ihit){
341  if(trkhits[ihit]->Plane() < plane){ continue; }
342  if(trkhits[ihit]->Plane() > plane){ break; }
343  // now have a hit on the same plane
344  if(trkhits[ihit]->Cell() == cellmax){ topok = true; }
345  else if(trkhits[ihit]->Cell() == cellmin){ botok = true; }
346  }
347  if(topok && botok){ tricells.push_back(allcellsontrack[icell]); }
348  }
349 
350  // loop over all the cells this track passes through
351  for(unsigned int icell = 0; icell < tricells.size(); ++icell){
352 
353  unsigned int cell = tricells[icell].Cell();
354  unsigned int plane = tricells[icell].Plane();
355 
356  // if plane is not less than plane max skip
357  if(plane > fplanemax-1){ continue; }
358 
359  // check if this channel is enabled
360  bool enable = true;
361  // need to double check this is working, leave disabled for now
362  // get dcm
363  // const daqchannelmap::lchan LogicalChannel = cmap->encodeLChan(fDet,plane,cell);
364  // // get the dcm id
365  // const uint32_t dcmid = cmap->getDCMOffline(LogicalChannel);
366  // // get the feb id
367  // const int febid = cmap->getDCMLinkOffline(LogicalChannel);
368  // // get the run history dcm object for this dcm id
369  // nova::dbi::RunHistory::DCM dcm = frunhist.GetDCM(dcmid);
370  // // get enable information for this febid
371  // enable = dcm.feb[febid].isEnabled;
372  // skip if not enabled
373  if(!enable){ continue; }
374 
375  // get the center of this cell
376  double xyz[3];
377  geom->Plane(plane)->Cell(cell)->GetCenter(xyz);
378 
379  double deltax,deltay,deltaz,w;
380  if(geom->Plane(plane)->View() == geo::kX){
381  deltax = 2.0*geom->Plane(plane)->Cell(cell)->HalfW();
382  deltay = (dy/dx)*deltax;
383  deltaz = (dz/dx)*deltax;
384  // find the w position (y) for this cell
385  w = (xyz[2]-iz)*dy/dz + iy;
386  }
387  else{
388  deltay = 2.0*geom->Plane(plane)->Cell(cell)->HalfW();
389  deltax = (dx/dy)*deltay;
390  deltaz = (dz/dy)*deltay;
391  // find the w position (x) for this cell
392  w = (xyz[2]-iz)*dx/dz + ix;
393  }
394 
395  double tripath = sqrt(deltax*deltax+deltay*deltay+deltaz*deltaz);
396  double apddist = geom->DistToElectronics(w,*(geom->Plane(plane)->Cell(cell)));
397 
398  // loop over hits and see if this is a hit or not
399  bool hit = false;
400  for(unsigned int ihit = 0; ihit < trkhits.size(); ++ihit){
401  if(trkhits[ihit]->Cell() == cell && trkhits[ihit]->Plane() == plane){
402  hit = true;
403  break;
404  }
405  }
406 
407  // Fill some histograms
408  int cellbin = 384*plane+cell;
409 
410  if(geom->Plane(plane)->View() == geo::kX){
411  cellsVwVlenx->Fill(w,tripath);
412  cellsVapddistVlenx->Fill(apddist,tripath);
413  }
414  else{
415  cellsVwVleny->Fill(w,tripath);
416  cellsVapddistVleny->Fill(apddist,tripath);
417  }
418  cellsVwcell->Fill(cellbin,w);
419  cellsVapddistcell->Fill(cellbin,apddist);
420 
421  if(hit){
422  hitsVwcell->Fill(cellbin,w);
423  hitsVapddistcell->Fill(cellbin,apddist);
424  if(geom->Plane(plane)->View() == geo::kX){
425  hitsVwVlenx->Fill(w,tripath);
426  hitsVapddistVlenx->Fill(apddist,tripath);
427  }
428  else{
429  hitsVwVleny->Fill(w,tripath);
430  hitsVapddistVleny->Fill(apddist,tripath);
431  }
432  }
433 
434  // make plots with path length cut
435  if(tripath > minpath && tripath < maxpath){
436  cellsVwcelllencut->Fill(cellbin,w);
437  cellsVapddistcelllencut->Fill(cellbin,apddist);
438  if(hit){
439  hitsVwcelllencut->Fill(cellbin,w);
440  hitsVapddistcelllencut->Fill(cellbin,apddist);
441  }
442  }
443 
444  } // end of loop over tricells
445 
446  } // end of loop over tracks
447 
448  } // end of analyze
449 
451  {
452  NEvents->SetBinContent(1,fnevents);
453  }
454 
456 
457 } // end of namespace
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
void analyze(art::Event const &e)
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
std::map< std::string, double > xmax
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
X or Y views.
Definition: PlaneGeo.h:30
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
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
double HalfW() const
Definition: CellGeo.cxx:191
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Definition: Run.h:31
void reconfigure(fhicl::ParameterSet const &p)
Double_t ymax
Definition: plot.C:25
void CountCellsOnLine(double X1, double Y1, double Z1, double X2, double Y2, double Z2, std::vector< OfflineChan > &Xhitsonline, std::vector< OfflineChan > &Yhitsonline)
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
CDPStorage service.
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
nova::dbi::RunHistory frunhist
size_type size() const
Definition: PtrVector.h:308
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
void beginRun(art::Run const &r)
T * make(ARGS...args) const
double DetHalfWidth() const
HitEfficiency(fhicl::ParameterSet const &p)
double DistToElectronics(double localz, const CellGeo &cell) const
get distance from local z position in cell to apd in cm, including pigtail
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
TRandom3 r(0)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
unsigned int NPlanes() const
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
Float_t track
Definition: plot.C:35
Float_t w
Definition: plot.C:20
int InstrumentedDetEnds(double &frontDwnstrm, double &backDwnstrm, double &frontUpstrm, double &backUpstrm)
give all detector ends information; most general
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124
enum BeamMode string