PEResponse_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: PEResponse
3 // Module Type: Analyzer
4 // File: PEResponse_module.cc
5 // \brief A module to study channel response at a very basic level.
6 //
7 ////////////////////////////////////////////////////////////////////////
8 
9 // art and FHiCL includes
19 #include "fhiclcpp/ParameterSet.h"
20 
21 // NOvASoft includes
22 #include "RecoBase/Track.h"
23 #include "Geometry/Geometry.h"
25 #include "RunHistory/RunHistory.h"
28 #include "DAQChannelMap/DAQChannelMap.h"
29 #include "NovaDAQConventions/DAQConventions.h"
31 
32 // C++
33 #include <memory>
34 #include <string>
35 #include <iostream>
36 #include <vector>
37 
38 // ROOT
39 #include "TFile.h"
40 #include "TTree.h"
41 #include "TVector3.h"
42 #include "TH1.h"
43 #include "TH2.h"
44 
45 ////////////////////////////////////////////////////////////////////////
46 ////////////////////////////////////////////////////////////////////////
47 
48 
49 namespace calib
50 {
51  class PEResponse : public art::EDAnalyzer
52  {
53  public:
54  explicit PEResponse ( fhicl::ParameterSet const& p );
55  virtual ~PEResponse ();
56 
57  void analyze ( art::Event const& e );
58  void reconfigure ( fhicl::ParameterSet const& p );
59  void beginJob ( );
60  void beginRun ( art::Run const& r );
61  void endSubRun ( art::SubRun const& sr);
62  private:
63 
64  void SetDefault();
66  std::string fTrkLabel; ///< Track module label
67  std::string fPCHitLabel; ///< PCHit module label
68  std::string fQualXYLabel; ///< Instance label for hits with xz quality
69  std::string fQualZLabel; ///< Instance label for hits with z quality
70  std::string fQualAvgLabel; ///< Instance label for hits with avg quality
71  std::string fQualBadLabel; ///< Instance label for hits with bad track quality
72  std::string fExpLabel; ///< Cosmic exposure module label name
73 
74  TTree* fHitTree; ///< Tree to save hit information
75  std::vector<TH2F*> fThreshHists;
76  TH1F* hx[15];
77  TH1F* hy[15];
78 
79  TH2F* hthresh;
80 
81  TH1F* hall;
82  TH1F* hexp;
83  //TH1F* ncellperplane;
84  float fPe;
85  float fW;
86  float fPath;
87  float fTns;
88  float fReadDist;
89  float fFlightDist;
90  bool fIsBad;
91  bool fIsXYAdj;
92  bool fIsZAdj;
93  bool fIsAvg;
94  bool fIsBadTrk;
95  //bool fHasApd;
96  //bool fIsEnabled;
97  int fEvent;
98  int fRun;
99  int fSubRun;
100  int fPlane;
101  int fCell;
102  int fDcm;
103  int fDiblock;
104  int fApd;
105  int fPixel;
106  int fView;
109 
111  nova::dbi::RunHistory fRunHist; //Run history object
112 
113 
114  }; // end of PEResponse class definition
115 
116  ////////////////////////////////////////////////////////////////////////
117  ////////////////////////////////////////////////////////////////////////
118 
119  ///< Constructor
121  :EDAnalyzer( pset )
122  {
123  this->reconfigure(pset);
124  }
125 
126  ////////////////////////////////////////////////////////////////////////
127  ////////////////////////////////////////////////////////////////////////
128 
129 
130  ///< Destructor
132  {
133 
134  }
135 
136  ////////////////////////////////////////////////////////////////////////
137  ////////////////////////////////////////////////////////////////////////
138 
139 
140  ///< Reconfigure parameter set
142  {
143  fTrkLabel = pset.get< std::string >("TrkLabel");
144  fExpLabel = pset.get< std::string >("ExpLabel");
145  fCellHitLabel = pset.get< std::string >("CellHitLabel");
146  fPCHitLabel = pset.get< std::string >("PCHitLabel");
147  fQualXYLabel = pset.get< std::string >("QualXYLabel");
148  fQualZLabel = pset.get< std::string >("QualZLabel");
149  fQualAvgLabel = pset.get< std::string >("QualAvgLabel");
150  fQualBadLabel = pset.get< std::string >("QualBadLabel");
151  fFillHitTree = pset.get< bool >("FillHitTree");
152  }
153 
154 
155  ////////////////////////////////////////////////////////////////////////
156  ////////////////////////////////////////////////////////////////////////
157 
158  ///< beginJob function
160  {
162 
163  // declare all the histogram necessary
164 
165  hall = tfs->make<TH1F>("hall", "All W; PE; Number of Hits", 1000, 0, 1000);
166  hexp = tfs->make<TH1F>("hexp", "Total Exposure;;",1,0,1);
167 
168 
169  for( int i = 0; i < 15; i++){
170 
171  char *hxname = new char[10];
172  sprintf(hxname, "hx%d", i);
173  char *hyname = new char[10];
174  sprintf(hyname, "hy%d", i);
175 
176  char *htitle = new char[100];
177  int offset = -750;
178  int lower = (offset+(i*100));
179  int upper = offset+(100+(i*100));
180  sprintf(htitle, "%d cm < W < %d cm ; PE; Number of Hits", lower, upper );
181 
182  hx[i] = tfs->make<TH1F>(hxname, htitle, 1000, 0, 1000);
183  hy[i] = tfs->make<TH1F>(hyname, htitle, 1000, 0, 1000);
184  }
185 
186 
187  fHitTree = tfs->make<TTree>("HitTree","Hit Tree");
188  // define branches of the tree
189  fHitTree->Branch("run", &fRun);
190  fHitTree->Branch("subrun", &fSubRun);
191  fHitTree->Branch("event", &fEvent);
192  fHitTree->Branch("plane", &fPlane);
193  fHitTree->Branch("cell", &fCell);
194  fHitTree->Branch("diblock", &fDiblock);
195  fHitTree->Branch("dcm", &fDcm);
196  fHitTree->Branch("apd", &fApd);
197  fHitTree->Branch("pixel", &fPixel);
198  fHitTree->Branch("view", &fView);
199  fHitTree->Branch("pe", &fPe);
200  fHitTree->Branch("w", &fW);
201  fHitTree->Branch("path", &fPath);
202  fHitTree->Branch("tns", &fTns);
203  fHitTree->Branch("readdist", &fReadDist);
204  fHitTree->Branch("flightdist", &fFlightDist);
205  fHitTree->Branch("isxyadj", &fIsXYAdj);
206  fHitTree->Branch("iszadj", &fIsZAdj);
207  fHitTree->Branch("isavg", &fIsAvg);
208  fHitTree->Branch("isbadtrk", &fIsBadTrk);
209  fHitTree->Branch("isbad", &fIsBad);
210 
211  }// end of beginJob
212 
213 
214  ////////////////////////////////////////////////////////////////////////
215  ////////////////////////////////////////////////////////////////////////
216  ///< beginRun function
218  {
220  int run = r.run();
223  nova::dbi::RunHistory runhist(run,det);
224 
225  fRunHist = runhist;
227 
228  char *htitle = new char[100];
229  sprintf(htitle, "Threshold Map For Run %d;Plane;Cell", run);
230 
231  char *hname = new char[20];
232  sprintf(hname, "thresh-%d", run);
233  hthresh = tfs->make<TH2F>(hname, htitle, 1792, -0.5, 1791.5, 384, -0.5, 383.5);
234 
235  int ndb = fRunHist.NDiBlocks();
236 
237  for( int idb = 0; idb < ndb; ++idb){
239  std::vector< nova::dbi::RunHistory::DCM > dcms = thisdb.dcm;
240  int ndcm = dcms.size();
241 
242  for( int idcm = 0; idcm < ndcm; ++idcm){
243  nova::dbi::RunHistory::DCM thisdcm = dcms[ idcm];
244  std::vector< nova::dbi::RunHistory::FEB > febs = thisdcm.feb;
245  int nfeb = febs.size();
246 
247  for(int ifeb =0; ifeb < nfeb; ++ifeb){
248  nova::dbi::RunHistory::FEB thisfeb = febs[ ifeb];
249  std::vector< int16_t > thresh = thisfeb.pixelThresh;
250  int npix = thresh.size();
251 
252  for(int ipix = 0; ipix < npix; ++ipix){
253 
254  uint32_t daqchan = cmap->encodeDChan( det, thisdb.num, idcm+1, ifeb, ipix);
255  uint32_t logchan = cmap->encodeLChan(daqchan);
256  int cell = cmap->getCell( logchan );
257  int plane = cmap->getPlane( logchan );
258  // If threshold is 4095, the pixel was masked out
259  // Don't fill the histogram for it, otherwise
260  // it dominates the colz scale.
261  if (thresh[ipix] != 4095)
262  hthresh->Fill( plane, cell, thresh[ipix]);
263 
264  }// end loop over pixels
265  }// end loop over febs
266  }// end loop over dcms
267  }// end loop over diblocks
268 
269  fThreshHists.push_back( hthresh );
270  }// end of beginRun
271 
272  ////////////////////////////////////////////////////////////////////////
273  ////////////////////////////////////////////////////////////////////////
274 
275  ///< analyze function
277  {
278  SetDefault();
279  fSubRun = evt.subRun();
280  fEvent = evt.event();
281  fRun = evt.run();
282 
283 
286 
288  evt.getByLabel( fTrkLabel, trkCol );
289 
290 
295 
296  int nTrks = trkCol->size();
297 
298 
299  for( int iTrk = 0; iTrk < nTrks; ++iTrk){
300 
301  art::Ptr< rb::Track > track(trkCol, iTrk);
302  // grab the pc hits associated with this track
303  if ((track->NXCell()==0) || (track->NYCell()==0)) continue;
304 
305  // get the pc hits associated with the cosmic track
306  std::vector<art::Ptr<caldp::PCHit> > pcxy = fmpcxy.at(iTrk);
307  std::vector<art::Ptr<caldp::PCHit> > pcz = fmpcz.at(iTrk);
308  std::vector<art::Ptr<caldp::PCHit> > pcav = fmpcav.at(iTrk);
309  std::vector<art::Ptr<caldp::PCHit> > pcbt = fmpcbt.at(iTrk);
310 
311  fTrkPlanes = track->ExtentPlane();
312  // Now we will loop over all the pchits of different types
313  fIsXYAdj = true;
314  fIsZAdj = false;
315  fIsAvg = false;
316  fIsBadTrk = false;
317 
318 
319  for( auto const& iPc: pcxy ){
320  fW = iPc->W();
321  fPath = iPc->Path();
322  fTns = iPc->TNS();
323  fReadDist = iPc->ReadDist();
324  fFlightDist= iPc->FlightLen();
325  fPlane = iPc->Plane();
326  fCell = iPc->Cell();
327  fDcm = iPc->DCM();
328  fDiblock = iPc->Diblock();
329  fApd = iPc->APD();
330  fPixel = iPc->Pixel();
331  fView = iPc->View();
332  fIsBad = badc->IsBad( fPlane, fCell);
333  fPe = iPc->PE();
334 
335  int offset = -750;
336  if( fPlane < 128 && fIsXYAdj && fTrkPlanes > 20){
337 
338  for( int iH = 0; iH < 15; ++iH){
339  hall->Fill(fPe);
340  if( fW > offset+(iH*100) && fW < offset+(100+(iH*100)) ){
341  if( fView== geo::kX )
342  hx[iH]->Fill( fPe);
343  if( fView== geo::kY )
344  hy[iH]->Fill(fPe);
345  }
346  }// end loop over iH
347  }
348 
349  fHitTree->Fill();
350  }// end loop over pcxy
351 
352 
353  fIsXYAdj = false;
354  fIsZAdj = true;
355  fIsAvg = false;
356  fIsBadTrk = false;
357  for( auto const& iPc: pcz ){
358  fW = iPc->W();
359  fPath = iPc->Path();
360  fTns = iPc->TNS();
361  fReadDist = iPc->ReadDist();
362  fFlightDist= iPc->FlightLen();
363  fPlane = iPc->Plane();
364  fCell = iPc->Cell();
365  fDcm = iPc->DCM();
366  fDiblock = iPc->Diblock();
367  fApd = iPc->APD();
368  fPixel = iPc->Pixel();
369  fView = iPc->View();
370  fIsBad = badc->IsBad( fPlane, fCell);
371  fPe = iPc->PE();
372  fHitTree->Fill();
373  }// end loop over pcz
374 
375 
376  fIsXYAdj = false;
377  fIsZAdj = false;
378  fIsAvg = true;
379  fIsBadTrk = false;
380  for( auto const& iPc: pcav ){
381  fW = iPc->W();
382  fPath = iPc->Path();
383  fTns = iPc->TNS();
384  fReadDist = iPc->ReadDist();
385  fFlightDist= iPc->FlightLen();
386  fPlane = iPc->Plane();
387  fCell = iPc->Cell();
388  fDcm = iPc->DCM();
389  fDiblock = iPc->Diblock();
390  fApd = iPc->APD();
391  fPixel = iPc->Pixel();
392  fView = iPc->View();
393  fIsBad = badc->IsBad( fPlane, fCell);
394  fPe = iPc->PE();
395  fHitTree->Fill();
396  }// end loop over pcav
397 
398 
399  fIsXYAdj = false;
400  fIsZAdj = false;
401  fIsAvg = false;
402  fIsBadTrk = true;
403  for( auto const& iPc: pcbt ){
404  fW = iPc->W();
405  fPath = iPc->Path();
406  fTns = iPc->TNS();
407  fReadDist = iPc->ReadDist();
408  fFlightDist= iPc->FlightLen();
409  fPlane = iPc->Plane();
410  fCell = iPc->Cell();
411  fDcm = iPc->DCM();
412  fDiblock = iPc->Diblock();
413  fApd = iPc->APD();
414  fPixel = iPc->Pixel();
415  fView = iPc->View();
416  fIsBad = badc->IsBad( fPlane, fCell);
417  fPe = iPc->PE();
418  fHitTree->Fill();
419  }// end loop over pcbt
420 
421  }// end loop over tracks
422 
423  }// end of analyze
424 
426 
427  fIsBad = false;
428  fIsXYAdj = false;
429  fIsZAdj = false;
430  fIsAvg = false;
431  fIsBadTrk= false;
432  fPe = -5;
433  fW = -5;
434  fPath = -5;
435  fTns = -5;
436  fReadDist= -5;
437  fEvent = -5;
438  fRun = -5;
439  fSubRun = -5;
440  fPlane = -5;
441  fCell = -5;
442  fDcm = -5;
443  fDiblock = -5;
444  fApd = -5;
445  fPixel = -5;
446  fView = -5;
447  fFlightDist= -5;
448  fTrkPlanes = -5;
449  }// end of SetDefault
450 
451 
453  {
455  sr.getByLabel(fExpLabel, exp);
456  if( exp.failedToGet())
457  return;
458  hexp->Fill( 0.5, exp->totlivetime );
459  }
460 
462 
463 }// end of namespace calib
SubRunNumber_t subRun() const
Definition: Event.h:72
int NDiBlocks()
gives number of active diblocks only, may be less than 14
Definition: RunHistory.h:394
void analyze(art::Event const &e)
< analyze function
std::string fExpLabel
Cosmic exposure module label name.
std::string fCellHitLabel
std::string fQualBadLabel
Instance label for hits with bad track quality.
bool LoadPixelInfo(int nAttempt=0)
Definition: RunHistory.cxx:903
const char * p
Definition: xmltok.h:285
DiBlock GetDiBlock(int i, bool loadAll=true)
get ith diblock is RH list (which only includes diblocks with activity), starting with i=0...
Vertical planes which measure X.
Definition: PlaneGeo.h:28
virtual ~PEResponse()
< Destructor
void beginRun(art::Run const &r)
DEFINE_ART_MODULE(TestTMapFile)
RunNumber_t run() const
Definition: Run.h:47
cell_t getCell(lchan logicalchan) const
Decode the cell number from an lchan.
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
std::vector< int16_t > pixelThresh
Definition: RunHistory.h:269
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Definition: Run.h:31
std::string fQualZLabel
Instance label for hits with z quality.
CDPStorage service.
std::vector< FEB > feb
Definition: RunHistory.h:291
static DAQChannelMap * getInstance(int detID)
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
int evt
std::string fPCHitLabel
PCHit module label.
nova::dbi::RunHistory fRunHist
std::string fTrkLabel
Track module label.
caf::StandardRecord * sr
void beginJob()
< beginJob function
EventNumber_t event() const
Definition: Event.h:67
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Definition: run.py:1
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
std::vector< DCM > dcm
Definition: RunHistory.h:311
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
void reconfigure(fhicl::ParameterSet const &p)
< Reconfigure parameter set
std::vector< TH2F * > fThreshHists
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
art::ServiceHandle< geo::Geometry > fGeom
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
plane_t getPlane(lchan logicalchan) const
Decode the plane number from an lchan.
std::string fQualAvgLabel
Instance label for hits with avg quality.
TRandom3 r(0)
cmap::CMap class source code
Definition: CMap.cxx:17
dchan encodeDChan(int detID, diblock_t diblock, dcm_id_t dcm, feb_t feb, pixel_t pixel) const
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
bool IsBad(int plane, int cell)
Float_t track
Definition: plot.C:35
void endSubRun(art::SubRun const &sr)
std::string fQualXYLabel
Instance label for hits with xz quality.
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
PEResponse(fhicl::ParameterSet const &p)
< Constructor
TTree * fHitTree
Tree to save hit information.