NuEEnergyAlg.cxx
Go to the documentation of this file.
1 
2 
7 
8 #include "RecoBase/Shower.h"
9 #include "RecoBase/CellHit.h"
10 #include "RecoBase/RecoHit.h"
11 #include "RecoBase/Vertex.h"
12 #include "ShowerLID/NuEEnergyAlg.h"
13 #include "Calibrator/Calibrator.h"
16 
17 #include "TMath.h"
18 #include <map>
19 
20 namespace slid{
21 
22  /*********************************************************************************************/
23  ////////////////////////////////////////////////////////////////////////
24 
26  {
27  this->reconfigure(pset);
28  mf::LogInfo("NuEEnergyAlg") << " NuEEnergyAlg::NuEEnergyAlg()\n";
29  fMapIsFilled = false;
30  fEvtId = art::EventID( 999,999,999);
31  }
32 
33  /*********************************************************************************************/
34  ////////////////////////////////////////////////////////////////////////
36  {
37 
38  }
39 
40  /*********************************************************************************************/
41  ////////////////////////////////////////////////////////////////////////
42 
44  {
45  fCorAbsCell = pset.get< bool >("CorAbsCell");
46  fUseStdCellE = pset.get< bool >("UseStdCellE");
47  fCorDataCell = pset.get< bool >("CorDataCell");
48  fUseUncalibHits = pset.get< bool >("UseUncalibHits");
49  fPEThreshold = pset.get< float >("PEThreshold");
50  }
51 
52  /*********************************************************************************************/
53  ////////////////////////////////////////////////////////////////////////
54 
55  double NuEEnergyAlg::CellEnergy( const rb::Shower * shower,
56  std::vector< const rb::Shower *> showercol,
57  const int & index,
58  const art::EventID & evtid,
59  bool useweight,
60  bool* isCalibrated
61  )
62 
63  {
65  const rb::CellHit * cHit = &(*(shower->Cell( index )));
66  if( fEvtId != evtid )
67  ClearVars( evtid );
68  if( isCalibrated )
69  *isCalibrated = true;
70  geo::View_t view = cHit->View();
71  double w = shower->W(cHit);
72 
73  //get distance to readout electronics
74  double xyzd[4];
75  calib::GetXYZD(cHit->OfflineChan(), w, xyzd);
76  double readoutD = xyzd[3];
77 
78  int detid = fGeom->DetId();
79 
80  rb::Prong prong = *(shower);
81 
82  rb::RecoHit rhit = shower->RecoHit( index );
83  double cellgev = 0;
84 
85  if(isCalibrated)
86  *isCalibrated = rhit.IsCalibrated();
87 
88  if(rhit.IsCalibrated()){
89  cellgev = rhit.GeV();
90  }
91  else if(fUseUncalibHits){
92 
93  double adc = double(shower->Cell( index )->ADC());
94  double Att = 1;
95  Att = CellADCToGeV( readoutD, view);
96  cellgev = adc/Att;
97  }
98  else
99  return cellgev;
100 
101 
102  double weight = shower->Weight( index );
103  if(useweight == false) weight = 1;
104 
105  double norm = 1;
106  if( cellgev !=0 )
107  norm = cHit->PE()/cellgev;
108 
109  GetHitShowerMap( showercol );
110  std::map< int, int > hitmap;
111  if( weight != 1 && useweight == true){
112  for( auto & cmap : fHitShowerMap){
113  if(cmap.first == shower->Cell( index )->Channel())
114  hitmap = cmap.second;
115  }// end loop over common cells
116 
117  norm = 0;
118 
119  for( auto & idxmap : hitmap ){
120  double w1 = showercol[ idxmap.first ]->W( showercol[ idxmap.first ]->Cell(idxmap.second).get() );
121 
122  rb::RecoHit rhit1 = cal->MakeRecoHit(*(showercol[ idxmap.first ]->Cell(idxmap.second)), w1);
123  double cellgev1 = 0;
124 
125  if(rhit1.IsCalibrated()){
126  cellgev1 = rhit1.GeV();
127  }
128  else if( fUseUncalibHits ){
129  art::Ptr<rb::CellHit> thishit = showercol[ idxmap.first ]->Cell(idxmap.second);
130  double adc1 = thishit->ADC();
131  double Att1 = 1;
132  double xyzd[4];
133  calib::GetXYZD(thishit->OfflineChan(), showercol[ idxmap.first ]->W(thishit.get()), xyzd);
134  double readoutD = xyzd[3];
135  Att1 = CellADCToGeV( readoutD, view);
136  cellgev1 = adc1/Att1;
137  }
138 
139  double ai = cellgev1/cHit->PE();
140  norm += showercol[ idxmap.first ]->Weight( idxmap.second )/ai;
141  }
142  }
143 
144  if( fUseStdCellE ){
145  if( rhit.IsCalibrated() )
146  return rhit.GeV()*weight/norm;
147  else{
148  if(isCalibrated)
149  *isCalibrated = false;
150  return 0;
151  }
152  }
153  return NonStdCellEnergy( cHit, detid, readoutD, w, norm, weight);
154 
155  }// end of CellEnergy
156 
157 
158  /*********************************************************************************************/
159  ////////////////////////////////////////////////////////////////////////
160 
162  const int & detid,
163  const double & readoutD,
164  const double & w,
165  const double & norm,
166  const double & weight)
167  {
169  geo::View_t view = cHit->View();
170  double gev = 0;
171 
172  if( detid != novadaq::cnv::kFARDET || !fCorAbsCell){
173  gev = cHit->PE() * weight/norm;
174  }
175 
176  else{
177  // for FD, Attenuation effect calibrated for development
178  // (Product S12.02.14), should not be changed to keep the
179  // energy scale identical to the one used in ANN training!
180 
181  int adc = cHit->ADC(0);
182 
183  double Att = 1;
184  if(view == geo::kX)
185  Att = CellADCToGeV( readoutD, view);
186  else if(view == geo::kY)
187  Att = CellADCToGeV( readoutD, view);
188  gev = adc*weight/Att;
189  }
190 
191  //gev = gev;
192  return gev;
193 
194  }// end of NonStdCellEnergy
195 
196 
197  /*********************************************************************************************/
198  ////////////////////////////////////////////////////////////////////////
199 
200  double NuEEnergyAlg::CellEDataMC( double cellEnergy, double w,
202 
203  {
204  double ecor = 1; //S12-11-16
205  int detid = fGeom->DetId();
206 
207  if( cellEnergy > 0 && w > 0 ){
208 
209  if( detid == novadaq::cnv::kFARDET ){
210 
211  if( view == geo::kX )
212  ecor = 0.707772+0.000573526*w;
213 
214  else if( view == geo::kY )
215  ecor = 0.707772+0.000573526*w;
216 
217  }
218 
219  else{
220 
221  if( view == geo::kX )
222  ecor = 0.707772+0.000573526*w;
223 
224  else if( view == geo::kY )
225  ecor = 0.707772+0.000573526*w;
226 
227  }
228  }
229 
230  if(ecor<0.0001)ecor = 1.0;
231  return ecor;
232 
233  }// end of CellEDataMC
234 
235  /*********************************************************************************************/
236  ////////////////////////////////////////////////////////////////////////
237 
238 
239  double NuEEnergyAlg::CellADCToGeV( double readoutD, geo::View_t view)
240  {
241  double Att = 4.08050e+03;
243  int detid = fGeom->DetId();
244 
245  if( detid == novadaq::cnv::kFARDET ){
246 
247  if( view==geo::kX )
248  Att = 4.08574e+03+3.49097e+04*TMath::Exp(-2.67297e-03*readoutD)
249  -1.62442e+04*TMath::Exp(-3.14363e-02*readoutD-(-3.93281e-05)
250  *readoutD*readoutD-1.49793e-08*readoutD*readoutD*readoutD);
251 
252  else if( view == geo::kY )
253  Att = 4.09347e+03+3.50384e+04*TMath::Exp(-2.68025e-03*readoutD)
254  -1.74280e+04*TMath::Exp(-3.08694e-02*readoutD-(-3.86870e-05)
255  *readoutD*readoutD-1.48737e-08*readoutD*readoutD*readoutD);
256  }
257 
258  else{
259 
260  if( view == geo::kX )
261 
262  Att = 4.08574e+03+3.49097e+04*TMath::Exp(-2.67297e-03*readoutD)
263  -1.62442e+04*TMath::Exp(-3.14363e-02*readoutD-(-3.93281e-05)
264  *readoutD*readoutD-1.49793e-08*readoutD*readoutD*readoutD);
265 
266  else if( view == geo::kY )
267  Att = 4.09347e+03+3.50384e+04*TMath::Exp(-2.68025e-03*readoutD)
268  -1.74280e+04*TMath::Exp(-3.08694e-02*readoutD-(-3.86870e-05)
269  *readoutD*readoutD-1.48737e-08*readoutD*readoutD*readoutD);
270 
271  }
272  return Att;
273 
274  }// end of CellADCToGeV
275 
276  /*********************************************************************************************/
277  ////////////////////////////////////////////////////////////////////////
278 
279  void NuEEnergyAlg::GetHitShowerMap( std::vector< const rb::Shower *> showercol)
280  {
281  if( MapFilled( showercol) )
282  return;
283 
284  //clear hit map in case it was used previously before filling
285  fHitShowerMap.clear();
286 
287  // loop over all the showers and find the cellhits common
288  // to 2 or more showers.
289  int nSh = showercol.size();
290 
291  for( int iSh = 0; iSh < nSh; ++iSh){
292 
293  art::PtrVector< rb::CellHit > ihits = showercol[iSh]->AllCells();
294  int nCells = ihits.size();
295  // loop over cells in ith shower
296  for( int iCell = 0; iCell < nCells; ++iCell ){
297 
298  // In the following loop, we will find which other
299  // showercol the ihits[iCell] belongs to.
300  for( int jSh = iSh+1; jSh < nSh; ++jSh){
301  art::PtrVector< rb::CellHit > jhits = showercol[jSh]->AllCells();
302  int mCells = jhits.size();
303  // loop over cells in jth shower
304  for(int jCell = 0 ; jCell < mCells; ++jCell){
305 
306  if( jhits[ jCell] == ihits[ iCell] ){
307  fHitShowerMap[ ihits[ iCell]->Channel() ][ iSh ] = iCell;
308  fHitShowerMap[ ihits[ iCell]->Channel() ][ jSh ] = jCell;
309  }
310  }// end loop over jshower cells
311  }// end j loop over showercol
312  }// end loop over ishower cells
313  }// end i loop over showercol
314  fShowerCol.clear();
315  fShowerCol = showercol;
316  fMapIsFilled = true;
317  return;
318  }
319 
320  /*********************************************************************************************/
321  ////////////////////////////////////////////////////////////////////////
322 
323  bool NuEEnergyAlg::MapFilled( std::vector< const rb::Shower *> showercol)
324  {
325 
326  if( !fMapIsFilled )
327  return fMapIsFilled;
328  // Check if this shower collection is the same as the one
329  // we already have a map for (may be).
330 
331  if( fShowerCol.size() != showercol.size() ){
332  fHitShowerMap.clear();
333  fShowerCol.clear();
334  fMapIsFilled = false;
335  return fMapIsFilled;
336  }
337 
338  for( int ish = 0, nsh = showercol.size(); ish < nsh ; ++ish){
339  if( !(fShowerCol[ish]->NCell() == showercol[ish]->NCell() &&
340  fShowerCol[ish]->Dir() == showercol[ish]->Dir() &&
341  fShowerCol[ish]->Start() == showercol[ish]->Start() )){
342  fHitShowerMap.clear();
343  fShowerCol.clear();
344  fMapIsFilled = false;
345  return fMapIsFilled;
346  }
347  }
348 
349  return fMapIsFilled;
350  }
351 
352  /*********************************************************************************************/
353  ////////////////////////////////////////////////////////////////////////
354 
356  std::vector< const rb::Shower *> showercol,
357  const art::Event& evt)
358  {
359  double showerDepE = 0;
360  int ncells = shower->NCell();
361  art::EventID evtid = evt.id();
362 
363  if( fEvtId != evtid )
364  ClearVars( evtid );
365 
366  for( int i = 0; i < ncells; ++i){
367 
368  double cellE = CellEnergy( shower, showercol, i, evtid);
369  showerDepE += cellE;
370  }
371 
372  return showerDepE;
373 
374  }// end of ShowerDepEnergy
375 
376 
377  /*********************************************************************************************/
378  ////////////////////////////////////////////////////////////////////////
379 
380  double NuEEnergyAlg::ECorrMC( double gev ){
381  //EM shower energy correction
382  double ecor = 5.68142e-01; // with pigtail
383 
384  if( fCorAbsCell == false )
385  ecor = 5.62343e-01; //no mc correction
386 
387  return ecor;
388  }// end of CorrMC
389 
390 
391  /*********************************************************************************************/
392  ////////////////////////////////////////////////////////////////////////
393 
394  double NuEEnergyAlg::ShowerEnergy( const rb::Shower * shower,
395  std::vector< const rb::Shower *> showercol,
396  const art::Event& evt)
397  {
398  double showerE = 0;
399  double showerDepE = ShowerDepEnergy( shower, showercol, evt);
400  double ecor = ECorrMC( showerDepE );
401  //TODO: Average position is calculated in the way of RecoJMShower. It is not clear whether it would be better
402  //to calculate the average position as commented out below. Results are slightly different, should be looked at.
403  //For now doing it this way for consistency.
404  double avgx = (shower->Start().X() + shower->Stop().X())/2.0;
405  double avgy = (shower->Start().Y() + shower->Stop().Y())/2.0;
406  //NOTE: Turn off energy weighting when calculating mean position. On cosmics the low PE hits can skew the
407  //average since attenuation is not taken into account.
408  //double avgx = shower->MeanX(rb::kByClusterWeightOnly);
409  //double avgy = shower->MeanY(rb::kByClusterWeightOnly);
410 
411  double poscor = 1.+6.28813e-05*avgx + 6.15165e-05*avgy; //with pigtail
412 
413  if( !fCorAbsCell )
414  poscor = 1.0+1.04866e-04*avgx + 6.40668e-05*avgy;
415 
416  showerE = showerDepE/(ecor*poscor);
417 
418  return showerE;
419 
420  }// end of ShowerEnergy
421 
422 
423  /*********************************************************************************************/
424  ////////////////////////////////////////////////////////////////////////
425 
426  double NuEEnergyAlg::HadronicDepEnergy( std::vector< const rb::Shower *> showercol,
427  const rb::Cluster * slice,
428  const int & elecShowerIndex,
429  const art::Event& evt)
430  {
432  if( fEvtId != evt.id() )
433  ClearVars( evt.id() );
434 
435  double hadDepE = 0;
436  int nShowers = showercol.size();
437  rb::Cluster excludedHits = *(slice);
438 
439  double meanx = slice->MeanX();
440  double meany = slice->MeanY();
441 
442  for( int i = 0; i < nShowers; ++i ){
443  const rb::Shower * shower = showercol[i];
444  int nhits = shower->NCell();
445 
446  for( int h = 0; h < nhits; ++h){
447  excludedHits.RemoveHit( shower->Cell(h) );
448  }
449  if ( i == elecShowerIndex )
450  continue;
451  else{
452  hadDepE += ShowerDepEnergy( showercol[i], showercol, evt );
453  }
454  }// end loop over showers
455 
456  int nExHits = excludedHits.NCell();
457  for( int i = 0; i < nExHits; ++i){
458  if( excludedHits.Cell(i)->PE() < fPEThreshold )
459  continue;
460  rb::RecoHit rHit;
461  if( excludedHits.Cell(i)->View() == geo::kX )
462  rHit = rb::RecoHit(cal->MakeRecoHit(*(excludedHits.Cell(i)), meany));
463  else
464  rHit = rb::RecoHit(cal->MakeRecoHit(*(excludedHits.Cell(i)), meanx));
465 
466  if( rHit.IsCalibrated() ){
467  hadDepE += rHit.GeV();
468  }
469  else if( fUseUncalibHits ){
470  art::Ptr< rb::CellHit > cHit = excludedHits.Cell(i);
471  double adc = double( cHit->ADC(0) );
472  double Att = 1;
473  geo::View_t view = fGeom->Plane( cHit->Plane() )->View();
474 
475  double xyzd[4];
476  if(view == geo::kX)
477  calib::GetXYZD(cHit->OfflineChan(), slice->MeanY(), xyzd);
478  else
479  calib::GetXYZD(cHit->OfflineChan(), slice->MeanX(), xyzd);
480 
481  Att = CellADCToGeV( xyzd[3], view);
482  hadDepE += adc/Att;
483  }
484 
485  }// end loop over orphan hits.
486  return hadDepE;
487 
488  }// end of HadronicDepEnergy
489 
490 
491  /*********************************************************************************************/
492  ////////////////////////////////////////////////////////////////////////
493 
494  double NuEEnergyAlg::HadronicEnergy( std::vector< const rb::Shower *> showercol,
495  const rb::Cluster * slice,
496  const int & elecShowerIndex,
497  const art::Event& evt)
498  {
499  double hadDepE = HadronicDepEnergy( showercol, slice, elecShowerIndex, evt);
500  double hadE = hadDepE/(4.00599e-01
501  - 5.08688e+00 * TMath::Exp(-1.07930e+01*hadDepE-(-1.17127e+01)*pow(hadDepE,2)-1.01452e+01*pow(hadDepE,3))
502  +4.71995e+00 * TMath::Exp(-1.02536e+01*hadDepE-(-9.81034e+00)*pow(hadDepE,2)-4.22648e+00*pow(hadDepE,3)));
503  return hadE;
504  }// end of HadronicEnergy
505 
506 
507  /*********************************************************************************************/
508  ////////////////////////////////////////////////////////////////////////
509 
510  double NuEEnergyAlg::NuEEnergy( std::vector< const rb::Shower *> showercol,
511  const rb::Cluster * slice,
512  const int & elecShowerIndex,
513  const art::Event& evt)
514  {
515  double elecE = ShowerEnergy( showercol[elecShowerIndex], showercol, evt );
516  double hadE = HadronicEnergy( showercol, slice, elecShowerIndex, evt);
517  double nuE = ( elecE + hadE );
518 
519  return nuE;
520 
521  }// end of NuEEnergy
522 
523  /*********************************************************************************************/
524  ////////////////////////////////////////////////////////////////////////
525 
526  double NuEEnergyAlg::VertexEnergy(const rb::Shower * shower,
527  std::vector< const rb::Shower *> showercol,
528  const rb::Vertex * vertex,
529  const rb::Cluster * slice,
530  const art::Event& evt)
531  {
534 
535  if( fEvtId != evt.id() )
536  ClearVars( evt.id() );
537 
538  TVector3 trkstart = shower->Start();
539  TVector3 trkstop = shower->Stop();
540 
541 
542  int startPlane, stopPlane, vtxPlane=0;
543 
545 
546  for(int trial = 0; trial < 4; ++trial){
547  // Figure out what cell the vertex was in. Step 3cm in z in case of
548  // failure. If that didn't work, it must be in plastic between cells,
549  // move laterally in x or y.
550  double offsetX = 0, offsetY = 0;
551  if(trial == 1) offsetX = -2;
552  if(trial == 2) offsetY = -2;
553  if(trial == 3) offsetX = offsetY = -2;
554 
555  cid = geom->CellId(vertex->GetX()+offsetX, vertex->GetY()+offsetY, vertex->GetZ(),
556  0, 0, 1, -1);
557  if(cid)
558  break;
559  }
560 
561  if( cid )
562  vtxPlane = geom->getPlaneID( cid );
563  else
564  vtxPlane = slice->MinPlane();
565 
566 
567  startPlane = shower->MinPlane();
568  stopPlane = shower->MaxPlane();
569 
570 
571  if( abs( stopPlane - vtxPlane) < abs( startPlane - vtxPlane)){
572  int plane = stopPlane;
573  stopPlane = startPlane;
574  startPlane = plane;
575  }
576  int maxPlane = std::max(stopPlane, startPlane);
577  int minPlane = std::min(stopPlane, startPlane);
578 
579  bool backtag = false;
580  if(stopPlane < startPlane)
581  backtag = true;
582 
583  rb::Cluster lesserSlice = (shower->NCell() < slice->NCell()) ? slice->Exclude( shower ) : rb::Cluster();
584 
585  art::PtrVector<rb::CellHit> sliceHits = lesserSlice.AllCells();
586 
587  int nsh = sliceHits.size();
588  int nShowers = showercol.size();
589 
590  for(int ic = 0; ic < nsh; ++ic){
591  bool matched = false;
592  for(int is = 0; is < nShowers; ++is){
593  int nhits = showercol[is]->NCell();
594  for( int ish = 0; ish < nhits; ++ish ){
595  if( *(sliceHits[ic]) == *(showercol[is]->Cell(ish)) ){
596  matched = true;
597  break;
598  }
599  }
600  }
601  // hits that do not belong to any showers are quite
602  // possibly noise hits and subjected to a PE cut.
603  if( !matched && sliceHits[ic]->PE() <= fPEThreshold )
604  lesserSlice.RemoveHit( sliceHits[ic]);
605  }
606 
607  sliceHits = lesserSlice.AllCells();
608  nsh = sliceHits.size();
609 
610  double vtxE = 0;
611 
612  for( int ish = 0; ish < nsh; ++ish){
613 
614  int dPlane = sliceHits[ish]->Plane() - startPlane;
615  if( backtag )
616  dPlane = -dPlane;
617  int plane = sliceHits[ish]->Plane();
618  int cell = sliceHits[ish]->Cell();
619 
620  if( ((dPlane >= -8 && dPlane < 0) ||
621  (dPlane < 8 && plane <= maxPlane && plane >= minPlane) ) ){
622 
623  geo::View_t view = sliceHits[ish]->View();
624  int detid = geom->DetId();
625 
626  double cellXYZ[3]={0,0,0};
627  geom->Plane(plane)->Cell(cell)->GetCenter(cellXYZ);
628  TVector3 trkdir = shower->Dir();
629 
630  double w;
631  if( view == geo::kY)
632  w = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
633  else
634  w = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
635 
636  if( fUseStdCellE ){
637  rb::RecoHit rhit(cal->MakeRecoHit(*(sliceHits[ish]), w));
638  if( rhit.IsCalibrated() )
639  vtxE += rhit.GeV();
640  }
641  else{
642  double readoutD;
643  double norm;
644 
645  readoutD = geom->DistToElectronics( w, *geom->Plane(plane)->Cell( cell ));
646  norm = sliceHits[ish]->PE()/cal->GetPECorr(*(sliceHits[ish]), w ) ;
647 
648  double cellgev = 0;
649  rb::RecoHit rhit = cal->MakeRecoHit( *(sliceHits[ish]), w);
650 
651  if(rhit.IsCalibrated()){
652  cellgev = rhit.GeV();
653  }
654  else if( fUseUncalibHits ){
655  double adc = double(sliceHits[ish]->ADC());
656  double Att = 1;
657  Att = CellADCToGeV( readoutD, view);
658  cellgev = adc/Att;
659  }
660 
661  norm = sliceHits[ish]->PE()/cellgev;
662 
663  double weight = 1;
664  double e = NonStdCellEnergy( sliceHits[ish].get(), detid, readoutD, w, norm, weight);
665  vtxE += e;
666  }
667  }
668  }// end loop over slice hits
669 
670  return vtxE;
671  }//End of VertexEnergy
672 
673  /*********************************************************************************************/
674  ////////////////////////////////////////////////////////////////////////
675 
676 
677  double NuEEnergyAlg::ShowerSumEnergy(std::vector< const rb::Shower *> showercol,
678  const art::Event& evt)
679  {
680  double sumE = 0;
681 
682  for( auto const iShower: showercol){
683  sumE += ShowerEnergy( iShower, showercol, evt);
684  }// end loop over showers
685  return sumE;
686 
687  }// End of ShowerSumEnergy
688 
689  /*********************************************************************************************/
690  ////////////////////////////////////////////////////////////////////////
691 
692 
693  void NuEEnergyAlg::ClearVars( const art::EventID& evtid ){
694 
695  fShowerCol.clear();
696  fHitShowerMap.clear();
697  fMapIsFilled = false;
698  fEvtId = evtid;
699 
700  }// End of ClearVars
701  /*********************************************************************************************/
702  ////////////////////////////////////////////////////////////////////////
703 
704 
705 }// end of jmshower namespace
double CellEDataMC(double gev, double w, geo::View_t view)
Correct cell DATA/MC difference in cell energy computation. Only used if fCorDataCell is true...
T max(const caf::Proxy< T > &a, T b)
std::vector< const rb::Shower * > fShowerCol
Definition: NuEEnergyAlg.h:144
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
double VertexEnergy(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const rb::Vertex *vertex, const rb::Cluster *slice, const art::Event &evt)
Returns energy within 8 planes of the start of the shower.
Definition: event.h:34
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
double NuEEnergy(const std::vector< const rb::Shower * > shower, const rb::Cluster *slice, const int &elecShowerIndex, const art::Event &evt)
Returns the neutrino energy, assuming &#39;shower&#39; is electron shower.
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Prong.cxx:135
void GetHitShowerMap(std::vector< const rb::Shower * > showercol)
Fills fHitShowerMap.
void reconfigure(const fhicl::ParameterSet &pset)
const Var weight
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
double GetX() const
Definition: Vertex.h:23
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
Definition: Cluster.h:47
void abs(TH1 *hist)
void RemoveHit(const art::Ptr< rb::CellHit > hit)
Remove hit from current cluster.
Definition: Cluster.cxx:290
double GetY() const
Definition: Vertex.h:24
const PlaneGeo * Plane(unsigned int i) const
geo::OfflineChan OfflineChan() const
Definition: CellHit.h:49
rb::Cluster Exclude(const rb::Cluster *excl) const
Create a cluster from this one, but with the hits of excl removed.
Definition: Cluster.cxx:233
virtual TVector3 Start() const
Definition: Prong.h:73
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
bool fCorDataCell
Should Data MC cell energy correction be used?
Definition: NuEEnergyAlg.h:126
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
double GetZ() const
Definition: Vertex.h:25
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
double MeanX(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:231
unsigned short Cell() const
Definition: CellHit.h:40
double ShowerSumEnergy(const std::vector< const rb::Shower * > shower, const art::Event &evt)
Returns sum of ShowerEnergy of all showers in showercol.
Far Detector at Ash River, MN.
double HadronicEnergy(const std::vector< const rb::Shower * > showercol, const rb::Cluster *slice, const int &elecShowerIndex, const art::Event &evt)
Returns the hadronic energy in the slice, assuming &#39;shower&#39; is the electron shower.
double MeanY(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:232
double NonStdCellEnergy(const rb::CellHit *cHit, const int &detid, const double &readoutD, const double &w, const double &norm, const double &weight)
float fPEThreshold
PE threshold for considering a hit in hadronic energy calculation.
Definition: NuEEnergyAlg.h:129
double HadronicDepEnergy(const std::vector< const rb::Shower * > showercol, const rb::Cluster *slice, const int &elecShowerIndex, const art::Event &evt)
Returns the deposited hadronic energy.
void ClearVars(const art::EventID &evtid)
Clears the privately owned variables of this class.
std::map< uint32_t, std::map< int, int > > fHitShowerMap
Definition: NuEEnergyAlg.h:138
double ShowerDepEnergy(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Returns the deposited energy of a shower. Does not contain corrections due to dead material or thresh...
T get(std::string const &key) const
Definition: ParameterSet.h:231
int getPlaneID(const CellUniqueId &id) const
double CellEnergy(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const int &index, const art::EventID &evtid, bool useweight=true, bool *isCalibrated=NULL)
Returns the energy deposited in a cell. Used in computation of total shower energy. Only used if fUseStdCellE is false.
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
TVector3 Stop() const
Endpoint of the shower.
Definition: Shower.cxx:55
int evt
float PE() const
Definition: CellHit.h:42
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
void GetXYZD(geo::OfflineChan chan, double w, double *xyzd)
Return position in world coordninates and distance to the readout.
Definition: CalibUtil.cxx:294
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
Vertex location in position and time.
Definition: View.py:1
size_type size() const
Definition: PtrVector.h:308
double ShowerEnergy(const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
Returns the total energy of a shower. along with corrections due to dead material and threshold effec...
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
bool fUseStdCellE
Should standard recoHit function, GeV(), be used to compute cell energy?
Definition: NuEEnergyAlg.h:128
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
float GeV() const
Definition: RecoHit.cxx:69
A Cluster with defined start position and direction.
Definition: Prong.h:19
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
Float_t norm
NuEEnergyAlg(fhicl::ParameterSet const &pset)
bool fUseUncalibHits
Do we want to use uncalibrated hits in the energy calculations?
Definition: NuEEnergyAlg.h:131
T const * get() const
Definition: Ptr.h:321
double DistToElectronics(double localz, const CellGeo &cell) const
get distance from local z position in cell to apd in cm, including pigtail
Calculates deposited and corrected energy of the electron shower and of electron flavoured neutrino...
A rb::Prong with a length.
Definition: Shower.h:18
art::EventID fEvtId
Event number.
Definition: NuEEnergyAlg.h:154
void geom(int which=0)
Definition: geom.C:163
double CellADCToGeV(double d, geo::View_t view)
Convert ADC to GeV based on JM attenuation fits.
uint32_t Channel() const
Definition: RawDigit.h:84
double ECorrMC(double gev)
EM shower energy correction for MC Only used if fCorAbsCell is false.
bool fMapIsFilled
Is fHitShowerMap filled for the current shower collection?
Definition: NuEEnergyAlg.h:151
int ncells
Definition: geom.C:124
cmap::CMap class source code
Definition: CMap.cxx:17
double GetPECorr(rb::CellHit const &cellhit, double w)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
Build slid::LID objects to store electron ID, if asked for, otherwise, calculate LID info and make av...
Definition: FillPIDs.h:13
art::ServiceHandle< geo::Geometry > fGeom
Definition: NuEEnergyAlg.h:124
bool MapFilled(std::vector< const rb::Shower * > showercol)
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20
EventID id() const
Definition: Event.h:56
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Definition: Cluster.cxx:209
bool fCorAbsCell
Should non-standard absolute energy calibration be used?
Definition: NuEEnergyAlg.h:127