CalHitAna_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 /// \brief A simple module to monitor the performance of Calhit.
3 /// \author Kanika Sachdev (ksachdev@physics.umn.edu)
4 /// \date Dec. 2012
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "RecoBase/CellHit.h"
9 #include "Geometry/Geometry.h"
11 
18 #include <fhiclcpp/ParameterSet.h>
21 
22 #include <vector>
23 #include <cmath>
24 
25 #include "TTree.h"
26 
27 namespace calhit{
28 
29  class CalHitAna : public art::EDAnalyzer {
30 
31  public:
32  explicit CalHitAna(fhicl::ParameterSet const& pset);
33  virtual ~CalHitAna();
34 
35  void analyze(art::Event const& evt);
36  void beginJob();
37  void reconfigure(const fhicl::ParameterSet& pset);
38 
39  private:
40 
42 
43  TTree* fHitTree; ///< Tree that gets filled for every hit
44  TTree* fNonNoiseHitTree;///< Tree that gets filled for every hit
45  TTree* fEvtTree; ///< Tree that gets filled for every event
46 
47  float fPe; ///< PE of the CellHit
48  int fAdc; ///< ADC0 of the CellHit
49  int fTdc; ///< TDC0 of the CellHit
50  float fTns; ///< TNS of the CellHit
51  float fFLSTime; ///< Energy weighted time of energy deposit in cell
52  float fReadDist; ///< Distance to readout APD from energy deposit (with pigtail)
53  float fFLSEnergy; ///< FLS Energy in this hit
54  int fIsNoise; ///< Is this hit a noise?
55  int fIsGoodTiming; ///< Good timing for this hit?
56  float fNFLSHits; ///< Number of FLSHits in this hit
57  unsigned int fPlane; ///< Plane of the CellHit
58  unsigned int fCell; ///< Cell of the CellHit
59 
60  // event level variables
61 
62  int fNoise; ///< Number of noise hits in event
63  int fPhysics; ///< Number of physics hits in event
64  int fAllHits; ///< Number of all hits in event
65  float fPePerTrueE; ///< PE per True Energy in event
66  float fAdcPerTrueE; ///< ADC per True Energy in event
67 
68  std::string fCellHitLabel; ///< Module label name from where to grab CellHits
69  std::string fSlicerLabel; ///< Module label for slices
70  };
71 }
72 
73 /*********************************************************************************************/
74 ////////////////////////////////////////////////////////////////////////
75 
76 namespace calhit {
77 
78 
80  : EDAnalyzer(pset)
81  {
82  this->reconfigure(pset);
83  }
84 
85 
86  /*********************************************************************************************/
87  ////////////////////////////////////////////////////////////////////////
88 
89 
91 
92 
93  /*********************************************************************************************/
94  ////////////////////////////////////////////////////////////////////////
95 
96 
98  {
99  fCellHitLabel = pset.get<std::string>("CellHitLabel");
100  fSlicerLabel = pset.get<std::string>("SlicerLabel");
101  }
102 
103  /*********************************************************************************************/
104  ////////////////////////////////////////////////////////////////////////
105 
106 
108 
110 
111  fHitTree = tfs->make<TTree>("fHitTree","fHitTree");
112 
113  fHitTree->Branch("Pe", &fPe);
114  fHitTree->Branch("Adc", &fAdc);
115  fHitTree->Branch("Tdc", &fTdc);
116  fHitTree->Branch("Tns" , &fTns);
117  fHitTree->Branch("FlsTime", &fFLSTime);
118  fHitTree->Branch("ReadDist", &fReadDist);
119  fHitTree->Branch("Plane", &fPlane);
120  fHitTree->Branch("Cell", &fCell);
121  fHitTree->Branch("NFlsHits", &fNFLSHits);
122  fHitTree->Branch("FlsEnergy", &fFLSEnergy);
123  fHitTree->Branch("IsNoise", &fIsNoise);
124  fHitTree->Branch("IsGoodTiming", &fIsGoodTiming);
125 
126  fNonNoiseHitTree = tfs->make<TTree>("fNonNoiseHitTree","fNonNoiseHitTree");
127 
128  fNonNoiseHitTree->Branch("Pe", &fPe);
129  fNonNoiseHitTree->Branch("Adc", &fAdc);
130  fNonNoiseHitTree->Branch("Tdc", &fTdc);
131  fNonNoiseHitTree->Branch("Tns" , &fTns);
132  fNonNoiseHitTree->Branch("FlsTime", &fFLSTime);
133  fNonNoiseHitTree->Branch("ReadDist", &fReadDist);
134  fNonNoiseHitTree->Branch("Plane", &fPlane);
135  fNonNoiseHitTree->Branch("Cell", &fCell);
136  fNonNoiseHitTree->Branch("NFlsHits", &fNFLSHits);
137  fNonNoiseHitTree->Branch("FlsEnergy", &fFLSEnergy);
138  fNonNoiseHitTree->Branch("IsNoise", &fIsNoise);
139  fNonNoiseHitTree->Branch("IsGoodTiming", &fIsGoodTiming);
140 
141  fEvtTree = tfs->make<TTree>("fEvtTree","fEvtTree");
142 
143  fEvtTree->Branch("Noise", &fNoise);
144  fEvtTree->Branch("PePerTrueE", &fPePerTrueE);
145  fEvtTree->Branch("AdcPerTrueE",&fAdcPerTrueE);
146  fEvtTree->Branch("Physics", &fPhysics);
147  fEvtTree->Branch("AllHits", &fAllHits);
148  }
149 
150  /*********************************************************************************************/
151  ////////////////////////////////////////////////////////////////////////
152 
153 
154 
156 
158 
160  evt.getByLabel(fCellHitLabel, cHitCol);
161 
163  evt.getByLabel(fSlicerLabel, sliceHandle);
164 
165  // Set default values
166 
167  fPlane = -5;
168  fCell = -5;
169  fPe = -5;
170  fTdc = -5;
171  fTns = -5;
172  fAdc = -5;
173  fNoise = -5;
174  fPePerTrueE = -5;
175 
177 
178  // Fill the hit tree
179  for(int i = 0; i < (int)cHitCol->size(); ++i){
180  art::Ptr< rb::CellHit> hit(cHitCol, i);
181  fPlane = hit->Plane();
182  fCell = hit->Cell();
183  fPe = hit->PE();
184  fAdc = hit->ADC();
185  fTdc = hit->TDC();
186  fTns = hit->TNS();
187  fIsGoodTiming = hit->GoodTiming();
188 
189  if(bt->HaveTruthInfo())
190  FillTruthInfo(hit);
191  else{
192  fFLSEnergy = -5;
193  fIsNoise = -5;
194  fFLSTime = -5;
195  fReadDist = -5;
196  fNFLSHits = -5;
197  }
198  fHitTree->Fill();
199  }
200 
201  // Fill tree for non-noise slices
202  int nsli = sliceHandle->size();
203  for(int isli = 0; isli < nsli; isli++){
204  if(sliceHandle->at(isli).IsNoise() )
205  continue;
206 
207  // Fill the non-noise hit tree
208  int nhits = sliceHandle->at(isli).NCell();
209  for(int ic = 0; ic < nhits; ++ic){
210  art::Ptr< rb::CellHit> hit = sliceHandle->at(isli).Cell(ic);
211  fPlane = hit->Plane();
212  fCell = hit->Cell();
213  fPe = hit->PE();
214  fAdc = hit->ADC();
215  fTdc = hit->TDC();
216  fTns = hit->TNS();
217  fIsGoodTiming = hit->GoodTiming();
218 
219  if(bt->HaveTruthInfo())
220  FillTruthInfo(hit);
221  else{
222  fFLSEnergy = -5;
223  fIsNoise = -5;
224  fFLSTime = -5;
225  fReadDist = -5;
226  fNFLSHits = -5;
227  }
228  fNonNoiseHitTree->Fill();
229  }// end loop over slice hits
230 
231  }// end loop over slices
232  // If truth information is available, fill event tree
233  if( bt->HaveTruthInfo()){
234  fNoise = 0;
235  fPePerTrueE = 0;
236  float totPe = 0;
237  float totAdc= 0;
238  float totE = 0;
239  fAllHits = cHitCol->size();
240 
241  for(int i = 0; i < fAllHits; ++i){
242  art::Ptr< rb::CellHit> hit(cHitCol, i);
243 
244  if( bt->IsNoise(hit) ){
245  fNoise += 1;
246  continue;
247  }
248  else
249  fPhysics += 1;
250 
251  const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(hit);
252  for(int j = 0; j < (int)flsHits.size(); j++){
253  totE += flsHits[j].GetEdep();
254  }//End loop over flshits in hit
255 
256  totPe += hit->PE();
257  totAdc+= hit->ADC();
258  }// End loop over hits
259 
260  if( totE != 0){
261  fPePerTrueE = totPe/totE;
262  fAdcPerTrueE = totAdc/totE;
263  }
264  fEvtTree->Fill();
265  }// EndIf HaveTruthInfo
266 
267  }// End of analyser
268 
269  //-------------------------------------------------------------------
270 
272  {
274  const std::vector< sim::FLSHit> flsHits = bt->HitToFLSHit(hit);
275  fNFLSHits = flsHits.size();
276  fIsNoise = bt->IsNoise(hit);
277  fFLSEnergy= 0;
278  fFLSTime = 0;
279 
280  // calculate the total FLS energy, energy weighted time
281  // and distance from readout.
282  double w = 0.0;
283  for( int ifls = 0; ifls < fNFLSHits; ifls++){
284  fFLSEnergy += flsHits[ifls].GetEdepBirks();
285  fFLSTime += flsHits[ifls].GetEdepBirks()*
286  (flsHits[0].GetEntryT() +flsHits[0].GetExitT())/2.0;
287 
288  if (hit->View() == geo::kX)
289  w = flsHits[ifls].GetEdepBirks()*
290  (flsHits[0].GetEntryY() + flsHits[0].GetExitY())/2.0;
291 
292  if (hit->View() == geo::kY)
293  w = flsHits[ifls].GetEdepBirks()*
294  (flsHits[0].GetEntryX() + flsHits[0].GetExitX())/2.0;
295 
296  }// end loop over FLS Hits
297 
298  if(fFLSEnergy > 0){
300  double xyzd[4];
301  w = w/fFLSEnergy;
302  calib::GetXYZD(hit->OfflineChan(),w,xyzd);
303  fReadDist = xyzd[3];
304 
305  }
306 
307  }
308 
309  //-------------------------------------------------------------------
310 
312 
313 } //End of namespace calhit
314 
315 ////////////////////////////////////////////////////////////////////////
float TNS() const
Definition: CellHit.h:46
unsigned int fPlane
Plane of the CellHit.
back track the reconstruction to the simulation
float fPe
PE of the CellHit.
int fIsNoise
Is this hit a noise?
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
int fTdc
TDC0 of the CellHit.
void reconfigure(const fhicl::ParameterSet &pset)
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
unsigned int fCell
Cell of the CellHit.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
TTree * fEvtTree
Tree that gets filled for every event.
float fReadDist
Distance to readout APD from energy deposit (with pigtail)
int fIsGoodTiming
Good timing for this hit?
std::string fCellHitLabel
Module label name from where to grab CellHits.
geo::OfflineChan OfflineChan() const
Definition: CellHit.h:49
DEFINE_ART_MODULE(TestTMapFile)
int fAdc
ADC0 of the CellHit.
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
unsigned short Cell() const
Definition: CellHit.h:40
float fFLSEnergy
FLS Energy in this hit.
TTree * fHitTree
Tree that gets filled for every hit.
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
T get(std::string const &key) const
Definition: ParameterSet.h:231
float PE() const
Definition: CellHit.h:42
void GetXYZD(geo::OfflineChan chan, double w, double *xyzd)
Return position in world coordninates and distance to the readout.
Definition: CalibUtil.cxx:294
int fNoise
Number of noise hits in event.
const double j
Definition: BetheBloch.cxx:29
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
void analyze(art::Event const &evt)
float fTns
TNS of the CellHit.
float fFLSTime
Energy weighted time of energy deposit in cell.
float fPePerTrueE
PE per True Energy in event.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
CalHitAna(fhicl::ParameterSet const &pset)
std::string fSlicerLabel
Module label for slices.
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: event.h:1
void geom(int which=0)
Definition: geom.C:163
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
This module creates CellHit object.
void FillTruthInfo(const art::Ptr< rb::CellHit > &hit)
float fNFLSHits
Number of FLSHits in this hit.
int fPhysics
Number of physics hits in event.
TTree * fNonNoiseHitTree
Tree that gets filled for every hit.
float fAdcPerTrueE
ADC per True Energy in event.
bool GoodTiming() const
Definition: CellHit.h:48
int fAllHits
Number of all hits in event.
Float_t w
Definition: plot.C:20
Encapsulate the geometry of one entire detector (near, far, ndos)