All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Types | Public Member Functions | Public Attributes | Private Attributes | List of all members
hit::ICARUSHitFinder Class Reference
Inheritance diagram for hit::ICARUSHitFinder:

Public Types

using ICARUSPeakFitParams_t = struct ICARUSPeakFitParams{float peakCenter
 
using ICARUSPeakParamsVec = std::vector< ICARUSPeakFitParams_t >
 

Public Member Functions

 ICARUSHitFinder (fhicl::ParameterSet const &pset)
 
void produce (art::Event &evt)
 
void beginJob ()
 
void endJob ()
 
void reconfigure (fhicl::ParameterSet const &p)
 
void expandHit (reco_tool::ICandidateHitFinder::HitCandidate &h, std::vector< float > holder, std::vector< reco_tool::ICandidateHitFinder::HitCandidate > how)
 
void computeBestLocalMean (std::vector< reco_tool::ICandidateHitFinder::HitCandidate > h, std::vector< float > holder, reco_tool::ICandidateHitFinder::MergeHitCandidateVec how, float &localmean)
 
void findMultiPeakParameters (const std::vector< float > &, const reco_tool::ICandidateHitFinder::HitCandidateVec &, ICARUSPeakParamsVec &, double &, int &, int) const
 
void findLongPeakParameters (const std::vector< float > &, const reco_tool::ICandidateHitFinder::HitCandidateVec &, ICARUSPeakParamsVec &, double &, int &, int) const
 
double ComputeChiSquare (TF1 func, TH1 *histo) const
 
double ComputeNullChiSquare (std::vector< float >) const
 
void setWire (int i)
 

Public Attributes

float peakCenterError
 
float peakSigma
 
float peakSigmaError
 
float peakAmplitude
 
float peakAmplitudeError
 
float peakTauLeft
 
float peakTauLeftError
 
float peakTauRight
 
float peakTauRightError
 
float peakFitWidth
 
float peakFitWidthError
 
float peakSlope
 
float peakSlopeError
 
float peakBaseline
 
float peakBaselineError
 

Private Attributes

std::vector< double > localmeans
 
unsigned int fDataSize
 
art::InputTag fDigitModuleLabel
 
std::string fSpillName
 
std::string fCalDataModuleLabel
 
std::string fHitLabelName
 
int fThetaAngle
 
bool fUncompressWithPed
 
std::unique_ptr
< reco_tool::ICandidateHitFinder
fHitFinderTool
 For finding candidate hits. More...
 
size_t fMaxMultiHit
 maximum hits for multi fit More...
 
double fChi2NDF
 
std::vector< int > fLongMaxHitsVec
 maximum Chisquared / NDF allowed for a hit to be saved More...
 
std::vector< int > fLongPulseWidthVec
 Sets width of hits used to describe long pulses. More...
 
TH1F * fFirstChi2
 
TH1F * fNullChi2
 
TH1F * fChi2
 
TH1F * fHeightC
 
TH1F * fWidthC
 
TH1F * fNoiseC
 
TH1F * fAreaC
 
TH1F * fnhwC
 
TH1F * fIntegralC
 
TH1F * fAreaInt
 
TH1F * fBaselineC
 
TH1F * fWire2771
 
TH1F * fHeightI2
 
TH1F * fWidthI2
 
TH1F * fNoiseI2
 
TH1F * fHeightI1
 
TH1F * fWidthI1
 
TH1F * fNoiseI1
 
double fMinWidth
 minimum initial width for ICARUS fit More...
 
double fMaxWidthMult
 multiplier for max width for ICARUS fit More...
 
int fFittingRange
 semi-width of interval where to fit hit More...
 
int fIntegratingRange
 semi-width of interval where to integrate fitting function More...
 
int iWire
 
ICARUShitFitCache fFitCache
 Cached functions for multi-peak fits. More...
 
ICARUSlongHitFitCache fLongFitCache
 Cached functions for long hits. More...
 
const geo::GeometryCorefGeometry = lar::providerFrom<geo::Geometry>()
 

Detailed Description

Definition at line 123 of file ICARUSHitFinder_module.cc.

Member Typedef Documentation

using hit::ICARUSHitFinder::ICARUSPeakFitParams_t = struct ICARUSPeakFitParams { float peakCenter

Definition at line 139 of file ICARUSHitFinder_module.cc.

Definition at line 156 of file ICARUSHitFinder_module.cc.

Constructor & Destructor Documentation

hit::ICARUSHitFinder::ICARUSHitFinder ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 237 of file ICARUSHitFinder_module.cc.

237  : EDProducer{pset}
238  {
239  this->reconfigure(pset);
240 
241  //LET HITCOLLECTIONCREATOR DECLARE THAT WE ARE GOING TO PRODUCE
242  //HITS AND ASSOCIATIONS TO RAW DIGITS BUT NOT ASSOCIATIONS TO WIRES
243  //(WITH NO PARTICULAR PRODUCT LABEL).
245  /*instance_name*/"");
246  }
void reconfigure(fhicl::ParameterSet const &p)
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:248

Member Function Documentation

void hit::ICARUSHitFinder::beginJob ( )

Definition at line 272 of file ICARUSHitFinder_module.cc.

273  {
274  // get access to the TFile service
275  art::ServiceHandle<art::TFileService> tfs;
276 
277 
278  // ======================================
279  // === Hit Information for Histograms ===
280  fFirstChi2 = tfs->make<TH1F>("fFirstChi2", "#chi^{2}", 10000, 0, 1000);
281  fNullChi2 = tfs->make<TH1F>("fNullChi2", "#chi^{2}", 100, 0, 10);
282 
283  fChi2 = tfs->make<TH1F>("fChi2", "#chi^{2}", 10000, 0, 1000);
284  fHeightC = tfs->make<TH1F>("fHeightC", "height(ADC#)", 100, 0, 100);
285  fWidthC = tfs->make<TH1F>("fWidthC", "width(samples)", 200, 0, 200);
286  fNoiseC = tfs->make<TH1F>("fNoiseC", "Noise Area(ADC#)", 100, 0, 100);
287  fAreaC = tfs->make<TH1F>("fAreaC", "", 150, 0, 1500);
288  fAreaC->GetXaxis()->SetTitle("Area(ADC#*ticks)");
289  fAreaC->GetYaxis()->SetTitle("events");
290 
291  fWire2771 = tfs->make<TH1F>("fWire2771", "", 4096, -0.5, 4095.5);
292  fWire2771->GetXaxis()->SetTitle("tick");
293  fWire2771->GetYaxis()->SetTitle("ADC#");
294 
295  fIntegralC = tfs->make<TH1F>("fIntegralC", "", 500, 0, 3000);
296  fIntegralC->GetXaxis()->SetTitle("Area(ADC#*ticks)");
297  fIntegralC->GetYaxis()->SetTitle("events");
298  fBaselineC = tfs->make<TH1F>("fBaselineC", "", 500, -10., 10.);
299  fBaselineC->GetXaxis()->SetTitle("Area(ADC#*ticks)");
300  fBaselineC->GetYaxis()->SetTitle("events");
301  fAreaInt = tfs->make<TH1F>("fAreaInt", "Area/Int", 100,0.,2.);
302 
303  fnhwC = tfs->make<TH1F>("fnhwC", "", 10, 0, 10);
304  fnhwC->GetXaxis()->SetTitle("n hits per wire");
305  fnhwC->GetYaxis()->SetTitle("wires");
306  fHeightI2 = tfs->make<TH1F>("fHeightI2", "height(ADC#)", 100, 0, 100);
307  fWidthI2 = tfs->make<TH1F>("fWidthI2", "width(samples)", 100, 0, 100);
308  fNoiseI2 = tfs->make<TH1F>("fNoiseI2", "Noise Area(ADC#)", 100, 0, 100);
309  fHeightI1 = tfs->make<TH1F>("fHeightI1", "height(ADC#)", 100, 0, 100);
310  fWidthI1 = tfs->make<TH1F>("fWidthI1", "width(samples)", 100, 0, 100);
311  fNoiseI1 = tfs->make<TH1F>("fNoiseI1", "Noise Area(ADC#)", 100, 0, 100);
312  // std::cout << " ICARUSHitfinder begin " << std::endl;
313  }
art::ServiceHandle< art::TFileService > tfs
void hit::ICARUSHitFinder::computeBestLocalMean ( std::vector< reco_tool::ICandidateHitFinder::HitCandidate h,
std::vector< float >  holder,
reco_tool::ICandidateHitFinder::MergeHitCandidateVec  how,
float &  localmean 
)

Definition at line 984 of file ICARUSHitFinder_module.cc.

985  {
986  const int bigw=130; //size of the window where to look for the minimum localmean value
987  const int meanw=70; //size of the window where the mean is calculated
988  const int outofbounds=9999;
989 
990  float samples1[bigw]; //list to contain samples bellow the startTick
991  float samples2[bigw]; //list to contain samples above the stopTick
992 
994 
995  float min1;
996  float min2;
997  int startTick=h.front().startTick;
998  int stopTick=h.back().stopTick;
999  float count1=0;
1000  float count2=0;
1001  unsigned int shift1=0;
1002  unsigned int shift2=0;
1003  int foundborder1=0;
1004  int foundborder2=0;
1005 
1006  // fill list of existing hits on this wire
1007  for(unsigned int j=0;j<how.size();j++)
1008  {
1009  std::vector<reco_tool::ICandidateHitFinder::HitCandidate> h2=how[j];
1010  if(h2.front().startTick!=h.front().startTick)
1011  hlist.push_back(h2);
1012  }
1013 
1014  // fill the arrays of samples to be examined
1015  for(unsigned int i=0;i<bigw;i++)
1016  {
1017  // remove samples from other hits in the wire
1018  for(unsigned int j=0;j<hlist.size();j++)
1019  {
1020  std::vector<reco_tool::ICandidateHitFinder::HitCandidate> h2=hlist[j];
1021  if(startTick-i-shift1 >= h2.front().startTick && startTick-i-shift1 <= h2.back().stopTick)
1022  shift1+=h2.back().stopTick-h2.front().startTick+1;
1023  else if(stopTick+i+shift2 >= h2.front().startTick && stopTick+i+shift2 <= h2.back().stopTick)
1024  shift2+=h2.back().stopTick-h2.front().startTick+1;
1025  }
1026 
1027  // fill the lists
1028  //if(startTick-i-shift1>=0)
1029  samples1[i]=holder[h.front().startTick-i-shift1];
1030  //else
1031  //samples1[i]=outofbounds;
1032 
1033  if(stopTick+i+shift2<=4095)
1034  samples2[i]=holder[h.back().stopTick+i+shift2];
1035  else
1036  samples2[i]=outofbounds;
1037  }
1038 
1039  // initialize counters
1040  for(int j=0;j<meanw;j++)
1041  {
1042  if(samples1[j]==outofbounds) foundborder1=1;
1043  if(samples2[j]==outofbounds) foundborder2=1;
1044 
1045  if(!foundborder1) count1+=samples1[j];
1046  if(!foundborder2) count2+=samples2[j];
1047  }
1048  min1=count1;
1049  min2=count2;
1050 
1051  //look for the local minima from the filled lists
1052  for(int j=0;j<bigw-meanw;j++)
1053  {
1054  if(count1<min1) min1=count1;
1055  if(samples1[j+meanw]==outofbounds) foundborder1=1;
1056  if(!foundborder1) count1+=samples1[j+meanw]-samples1[j];
1057 
1058  if(count2<min2) min2=count2;
1059  if(samples2[j+meanw]==outofbounds) foundborder2=1;
1060  if(!foundborder2) count2+=samples2[j+meanw]-samples2[j];
1061 
1062  if(foundborder1 && foundborder2) break;
1063  }
1064  if(count1<min1) min1=count1;
1065  if(count2<min2) min2=count2;
1066 
1067  // std::cout << " localmean count1 " << count1 << " count2 " << count2 << std::endl;
1068  //for the moment take the highest value
1069  if((foundborder1 && !foundborder2) || (shift1 && !shift2))
1070  localmean=((float) min2)/meanw;
1071  else if((!foundborder1 && foundborder2) || (!shift1 && shift2))
1072  localmean=((float) min1)/meanw;
1073  else
1074  localmean=(min1>min2)?((float) min1)/meanw :((float) min2)/meanw;
1075 
1076  //for the moment take the average value
1077  // h->localmean=(min1+min2)/(2.*meanw);
1078  }
std::vector< HitCandidateVec > MergeHitCandidateVec
double hit::ICARUSHitFinder::ComputeChiSquare ( TF1  func,
TH1 *  histo 
) const

Definition at line 1390 of file ICARUSHitFinder_module.cc.

1391  {
1392  double chi=0;
1393  int nb=histo->GetNbinsX();
1394  double wb=histo->GetBinWidth(0);
1395 
1396  int jp;
1397  for( jp=1;jp<nb;jp++) {
1398  if(histo->GetBinContent(jp)==0) break;
1399  double xb=jp*wb;
1400  double fv=(func)(xb);
1401  double hv=histo->GetBinContent(jp);
1402  double dv=hv-fv;
1403  double cv=dv/(histo->GetBinError(jp));
1404  chi+=cv*cv;
1405  //std::cout << " chi " << chi << std::endl;
1406 
1407  }
1408 // std::cout << " chi2mio" << chi << std::endl;
1409  //std::cout << " ndf " << ndf << std::endl;
1410  return chi/(jp-5);
1411  }
double hit::ICARUSHitFinder::ComputeNullChiSquare ( std::vector< float >  holder) const

Definition at line 1412 of file ICARUSHitFinder_module.cc.

1413  {
1414  double chi=0;
1415  int nb=33;
1416 
1417  int jp;
1418  for( jp=1;jp<nb;jp++) {
1419 
1420  double hv=holder[jp];
1421  double dv=hv;
1422  double cv=dv/2.4;
1423  chi+=cv*cv;
1424  // std::cout << " chi " << chi << std::endl;
1425 
1426  }
1427  // std::cout << " chi2null" << chi << std::endl;
1428  //std::cout << " ndf " << ndf << std::endl;
1429  return chi/(jp);
1430  }
void hit::ICARUSHitFinder::endJob ( )

Definition at line 315 of file ICARUSHitFinder_module.cc.

316  {
317  // std::cout << " ICARUSHitFinder endjob " << std::endl;
318 
319  }
void hit::ICARUSHitFinder::expandHit ( reco_tool::ICandidateHitFinder::HitCandidate h,
std::vector< float >  holder,
std::vector< reco_tool::ICandidateHitFinder::HitCandidate how 
)

Definition at line 896 of file ICARUSHitFinder_module.cc.

897  {
898  // Given a hit or hit candidate <hit> expand its limits to the closest minima
899  int nsamp=50;
900  int cut=1;
901  int upordown;
902  int found=0;
903 
904  unsigned int first=h.startTick;
905  unsigned int last =h.stopTick;
906 
907  std::vector<reco_tool::ICandidateHitFinder::HitCandidate>::iterator hiter;
908  std::vector<reco_tool::ICandidateHitFinder::HitCandidate> hlist;
909 
910  // fill list of existing hits on this wire
911  for(unsigned int j=0;j<how.size();j++)
912  {
914  if(h2.hitCenter!=h.hitCenter)
915  hlist.push_back(h2);
916  }
917 
918 
919  // look for first sample
920  while(!found)
921  {
922  if(first==0)
923  break;
924 
925  for(hiter=hlist.begin();hiter!=hlist.end();hiter++)
926  {
928  if(first==h2.stopTick)
929  found=1;
930  }
931 
932  if(found==1) break;
933 
934  upordown=0;
935  for(int l=0;l<nsamp/2;l++)
936  {
937  if(first-nsamp/2+l<4095) {
938  if(holder[first-nsamp/2+l+1]-holder[first-nsamp/2+l]>0) upordown++;
939  else if(holder[first-nsamp/2+l+1]-holder[first-nsamp/2+l]<0) upordown--;
940  } }
941  //std::cout << " checking " << first << " upordown " << upordown << std::endl;
942  if(upordown>cut)
943  first--;
944  else
945  found=1;
946  }
947 
948  // look for last sample
949  found=0;
950  while(!found)
951  {
952  if(last==4095)
953  {
954  found=1;
955  break;
956  }
957 
958  for(hiter=hlist.begin();hiter!=hlist.end();hiter++)
959  {
961  if(last==h2.startTick)
962  found=1;
963  }
964 
965  if(found==1) break;
966 
967  upordown=0;
968  for(int l=0;l<nsamp/2;l++)
969  {
970  if(last+nsamp/2-l>0) {
971  if(holder[last+nsamp/2-l]-holder[last+nsamp/2-l-1]>0) upordown++;
972  else if(holder[last+nsamp/2-l]-holder[last+nsamp/2-l-1]<0) upordown--;
973  } }
974 
975  if(upordown<-cut)
976  last++;
977  else
978  found=1;
979  }
980 
981  h.startTick=first;
982  h.stopTick=last;
983  }
void hit::ICARUSHitFinder::findLongPeakParameters ( const std::vector< float > &  roiSignalVec,
const reco_tool::ICandidateHitFinder::HitCandidateVec hitCandidateVec,
ICARUSPeakParamsVec peakParamsVec,
double &  chi2PerNDF,
int &  NDF,
int  iWire 
) const

Definition at line 1230 of file ICARUSHitFinder_module.cc.

1235  {
1236  TH1F* fHistogram=new TH1F("","",roiSignalVec.size(),0.,roiSignalVec.size());;
1237  std::string wireName = "PeakFitterHitSignal_" + std::to_string(iWire);
1238  fHistogram->SetName(wireName.c_str());
1239 
1240  if (hitCandidateVec.empty()) return;
1241 
1242  // in case of a fit failure, set the chi-square to infinity
1243  chi2PerNDF = std::numeric_limits<double>::infinity();
1244 
1245  int startTime = hitCandidateVec.front().startTick-fFittingRange;
1246  int endTime = hitCandidateVec.back().stopTick+fFittingRange;
1247  if(startTime<0) startTime=0;
1248  if(endTime>4095) endTime=4095;
1249 
1250  int roiSize = endTime - startTime;
1251 
1252  //std::cout << " roisize " << roiSize << std::endl;
1253 
1254  // Check to see if we need a bigger histogram for fitting
1255  if (roiSize > fHistogram->GetNbinsX())
1256  {
1257  std::string histName = "PeakFitterHitSignal_" + std::to_string(iWire);
1258  fHistogram = new TH1F(histName.c_str(),"",roiSize,0.,roiSize);
1259  fHistogram->Sumw2();
1260  }
1261 
1262  fHistogram->Reset();
1263  for(int idx = 0; idx < roiSize; idx++)
1264  fHistogram->SetBinContent(idx+1,roiSignalVec.at(startTime+idx));
1265 
1266  // Build the string to describe the fit formula
1267  std::string equation = "gaus(0)";
1268 
1269 
1270 
1271  // Now define the complete function to fit
1272  // TF1 Func("ICARUSfunc",fitlong,0,roiSize,1+7*hitCandidateVec.size());
1273  TF1& Func = *(fLongFitCache.Get(hitCandidateVec.size()));
1274  assert(&Func);
1275 
1276  // ### Setting the parameters for the ICARUS Fit ###
1277  Func.FixParameter(0,hitCandidateVec.size());
1278 
1279  int parIdx { 0 };
1280  for(auto const& candidateHit : hitCandidateVec)
1281  {
1282  double const peakMean = candidateHit.hitCenter - float(startTime);
1283  double const peakWidth = candidateHit.hitSigma;
1284  double const amplitude = candidateHit.hitHeight;
1285 
1286  Func.SetParameter(1+parIdx,0);
1287  Func.SetParameter(2+parIdx, amplitude);
1288  Func.SetParameter(3+parIdx, peakMean);
1289  Func.SetParameter(4+parIdx,peakWidth);
1290  Func.SetParameter(5+parIdx,peakWidth);
1291  Func.SetParameter(6+parIdx,2*peakWidth);
1292  Func.SetParameter(7+parIdx,0);
1293 
1294 
1295  Func.SetParLimits(1+parIdx, -5, 5);
1296  Func.SetParLimits(2+parIdx, 0.1 * amplitude, 10. * amplitude);
1297  Func.SetParLimits(3+parIdx, peakMean-peakWidth,peakMean+peakWidth);
1298  Func.SetParLimits(4+parIdx, std::max(fMinWidth, 0.01 * peakWidth), fMaxWidthMult * peakWidth);
1299  Func.SetParLimits(5+parIdx, std::max(fMinWidth, 0.01 * peakWidth), 4 * peakWidth);
1300  Func.SetParLimits(6+parIdx, 0,4*peakWidth);
1301  Func.SetParLimits(7+parIdx, -1,1);
1302 
1303  parIdx += 7;
1304 
1305  }
1306  int fitResult { -1 };
1307  try
1308  { fitResult = fHistogram->Fit(&Func,"QNWB","", 0., roiSize);
1309  }
1310  catch(...)
1311  {mf::LogWarning("GausHitFinder") << "Fitter failed finding a hit";}
1312 
1313  if(fitResult < -1)
1314  std::cout << " long fit cannot converge " << iWire << std::endl;
1315  // ##################################################
1316  // ### Getting the fitted parameters from the fit ###
1317  // ##################################################
1318  NDF = roiSize-7*hitCandidateVec.size();
1319  chi2PerNDF = (Func.GetChisquare() / NDF);
1320 
1321  parIdx = 0;
1322  peakParamsVec.clear();
1323  for(size_t idx = 0; idx < hitCandidateVec.size(); idx++)
1324  {
1325  ICARUSPeakFitParams_t peakParams;
1326 
1327  peakParams.peakAmplitude = Func.GetParameter(2+parIdx);
1328  peakParams.peakAmplitudeError = Func.GetParError(2+parIdx);
1329  peakParams.peakCenter = Func.GetParameter(3+parIdx) + float(startTime);
1330  peakParams.peakCenterError = Func.GetParError(3+parIdx);
1331 
1332  peakParams.peakTauRight = Func.GetParameter(4+parIdx);
1333  peakParams.peakTauRightError = Func.GetParError(4+parIdx);
1334  peakParams.peakTauLeft = Func.GetParameter(5+parIdx);
1335  peakParams.peakTauLeftError = Func.GetParError(5+parIdx);
1336  peakParams.peakFitWidth = Func.GetParameter(6+parIdx);
1337  peakParams.peakFitWidthError = Func.GetParError(6+parIdx);
1338  peakParams.peakSlope = Func.GetParameter(7+parIdx);
1339  peakParams.peakSlopeError = Func.GetParError(7+parIdx);
1340  peakParams.peakBaseline = Func.GetParameter(1+parIdx);
1341  peakParams.peakBaselineError = Func.GetParError(1+parIdx);
1342  // std::cout << " before adding peakparams size " << peakParamsVec.size() << std::endl;
1343  peakParamsVec.emplace_back(peakParams);
1344  // std::cout << " after adding peakparams size " << peakParamsVec.size() << std::endl;
1345 
1346  parIdx += 7;
1347 
1348  }
1349 
1350  //Func.Delete();
1351  bool writeWaveform=false;
1352  if(writeWaveform) {
1353  TFile *f = new TFile("fitICARUS.root","UPDATE");
1354  fHistogram->Write();
1355  f->Close();
1356  f->Delete();
1357  }
1358 
1359  fHistogram->Delete();
1360  return;
1361  }
struct ICARUSPeakFitParams{float peakCenter ICARUSPeakFitParams_t
double fMinWidth
minimum initial width for ICARUS fit
double fMaxWidthMult
multiplier for max width for ICARUS fit
ICARUSlongHitFitCache fLongFitCache
Cached functions for long hits.
int fFittingRange
semi-width of interval where to fit hit
virtual TF1 * Get(size_t nFunc)
Returns a function sum of nFunc base functions.
std::string to_string(WindowPattern const &pattern)
BEGIN_PROLOG could also be cout
void hit::ICARUSHitFinder::findMultiPeakParameters ( const std::vector< float > &  roiSignalVec,
const reco_tool::ICandidateHitFinder::HitCandidateVec hitCandidateVec,
ICARUSPeakParamsVec peakParamsVec,
double &  chi2PerNDF,
int &  NDF,
int  iWire 
) const

Definition at line 1079 of file ICARUSHitFinder_module.cc.

1084  {
1085 
1086  ICARUSPeakParamsVec peakParamsVec0;
1087  TH1F* fHistogram=new TH1F("","",roiSignalVec.size(),0.,roiSignalVec.size());;
1088 
1089  std::string wireName = "PeakFitterHitSignal_" + std::to_string(iWire);
1090  fHistogram->SetName(wireName.c_str());
1091 
1092  if (hitCandidateVec.empty()) return;
1093 
1094  // in case of a fit failure, set the chi-square to infinity
1095  chi2PerNDF = std::numeric_limits<double>::infinity();
1096 
1097  int startTime = hitCandidateVec.front().startTick-fFittingRange;
1098  int endTime = hitCandidateVec.back().stopTick+fFittingRange;
1099  if(startTime<0) startTime=0;
1100  if(endTime>4095) endTime=4095;
1101  int roiSize = endTime - startTime;
1102 
1103  //std::cout << " roisize " << roiSize << std::endl;
1104 
1105  // Check to see if we need a bigger hianchstogram for fitting
1106  if (roiSize > fHistogram->GetNbinsX())
1107  {
1108  std::string histName = "PeakFitterHitSignal_" + std::to_string(iWire);
1109  fHistogram = new TH1F(histName.c_str(),"",roiSize,0.,roiSize);
1110  //fHistogram->Sumw2();
1111  }
1112 
1113  fHistogram->Reset();
1114  for(int idx = 0; idx < roiSize; idx++)
1115  fHistogram->SetBinContent(idx+1,roiSignalVec.at(startTime+idx));
1116  for(int idx = 0; idx < roiSize; idx++)
1117  fHistogram->SetBinError(idx+1,2.4);
1118  // for(int idx = 0; idx < roiSize; idx++)
1119  // std::cout << " bin " << idx << " Error " << fHistogram->GetBinError(idx+1) << std::endl;
1120 
1121 
1122  // Now define the complete function to fit
1123  // TF1 Func("ICARUSfunc",fitf,0,roiSize,1+5*hitCandidateVec.size());
1124  TF1& Func = *(fFitCache.Get(hitCandidateVec.size()));
1125  assert(&Func);
1126 
1127  // ### Setting the parameters for the ICARUS Fit ###
1128  Func.FixParameter(0, hitCandidateVec.size());
1129 
1130  int parIdx{0};
1131  for(auto const& candidateHit : hitCandidateVec)
1132  {
1133  double const peakMean = candidateHit.hitCenter - float(startTime);
1134  double const peakWidth = candidateHit.hitSigma;
1135  // std::cout << " peakWidth " << peakWidth << std::endl;
1136  // std::cout << " hitcenter " << candidateHit.hitCenter << " starttime " << startTime <<std::endl;
1137 
1138 
1139  double const amplitude = candidateHit.hitHeight;
1140  // double meanLowLim = std::max(peakMean - fPeakRange * peakWidth, 0.);
1141  // double meanHiLim = std::min(peakMean + fPeakRange * peakWidth, double(roiSize));
1142 
1143  Func.SetParameter(1+parIdx,0);
1144  Func.SetParameter(2+parIdx, amplitude);
1145  Func.SetParameter(3+parIdx, peakMean);
1146  Func.SetParameter(4+parIdx,peakWidth);
1147  Func.SetParameter(5+parIdx,peakWidth);
1148 
1149  Func.SetParLimits(1+parIdx, -5, 5);
1150  Func.SetParLimits(2+parIdx, 0.1 * amplitude, 10. * amplitude);
1151  Func.SetParLimits(3+parIdx, peakMean-peakWidth,peakMean+peakWidth);
1152  Func.SetParLimits(4+parIdx, std::max(fMinWidth, 0.01 * peakWidth), fMaxWidthMult * peakWidth);
1153  Func.SetParLimits(5+parIdx, std::max(fMinWidth, 0.01 * peakWidth), fMaxWidthMult * peakWidth);
1154 
1155  parIdx += 5;
1156 
1157  }
1158 
1159  int fitResult(-1);
1160  // if(hitCandidateVec.size()>2) return;
1161  try
1162  { fitResult = fHistogram->Fit(&Func,"QNWB","", 0., roiSize);
1163  }
1164  catch(...)
1165  {mf::LogWarning("GausHitFinder") << "Fitter failed finding a hit";}
1166 
1167 
1168  // if(fitResult==0)
1169  // std::cout << " icarus fit converges " << iWire << std::endl;
1170  if(fitResult!=0)
1171 // std::cout << " icarus fit cannot converge " << iWire << std::endl;
1172  // ##################################################
1173  // ### Getting the fitted parameters from the fit ###
1174  // ##################################################
1175  NDF = roiSize-5*hitCandidateVec.size();
1176  chi2PerNDF = (Func.GetChisquare() / NDF);
1177 
1178  double chi2mio=ComputeChiSquare(Func,fHistogram);
1179 // std::cout << " chi2mio " << chi2mio << std::endl;
1180  chi2PerNDF=chi2mio;
1181  // for(int idx = 0; idx < roiSize; idx++)
1182  // std::cout << " bin " << idx << " Error " << fHistogram->GetBinError(idx+1) << std::endl;
1183 
1184 // std::cout << " chi2 " << Func.GetChisquare() << std::endl;
1185 // std::cout << " ndf " << NDF << std::endl;
1186 
1187  // std::cout << " chi2ndf " << chi2PerNDF<< std::endl;
1188  parIdx = 0;
1189  peakParamsVec0.clear();
1190 
1191  for(size_t idx = 0; idx < hitCandidateVec.size(); idx++)
1192  {
1193  ICARUSPeakFitParams_t peakParams;
1194 
1195  peakParams.peakAmplitude = Func.GetParameter(2+parIdx);
1196  peakParams.peakAmplitudeError = Func.GetParError(2+parIdx);
1197  peakParams.peakCenter = Func.GetParameter(3+parIdx) + float(startTime);
1198  peakParams.peakCenterError = Func.GetParError(3+parIdx);
1199  //std::cout << " rising time " << Func.GetParameter(3) << " falling time " <<Func.GetParameter(4) << std::endl;
1200  peakParams.peakTauRight = Func.GetParameter(4+parIdx);
1201  peakParams.peakTauRightError = Func.GetParError(4+parIdx);
1202  peakParams.peakTauLeft = Func.GetParameter(5+parIdx);
1203  peakParams.peakTauLeftError = Func.GetParError(5+parIdx);
1204  peakParams.peakBaseline = Func.GetParameter(1+parIdx);
1205  peakParams.peakBaselineError = Func.GetParError(1+parIdx);
1206  peakParams.peakFitWidth =0;
1207  peakParams.peakFitWidthError = 0;
1208  peakParams.peakSlope = 0;
1209  peakParams.peakSlopeError = 0;
1210  peakParamsVec.emplace_back(peakParams);
1211  parIdx += 5;
1212  //std::cout << " first center " << Func.GetParameter(2) + float(startTime) << std::endl;
1213  // std::cout << " second center " << Func.GetParameter(7) + float(startTime) << std::endl;
1214  }
1215 
1216  //Gaus.Delete();
1217  //Func.Delete();
1218  bool writeWaveform=false;
1219  if(writeWaveform) {
1220  TFile *f = new TFile("fitICARUS.root","UPDATE");
1221  fHistogram->Write();
1222  f->Close();
1223  f->Delete();
1224  }
1225 
1226  fHistogram->Delete();
1227  return;
1228  }
struct ICARUSPeakFitParams{float peakCenter ICARUSPeakFitParams_t
double fMinWidth
minimum initial width for ICARUS fit
double fMaxWidthMult
multiplier for max width for ICARUS fit
int fFittingRange
semi-width of interval where to fit hit
virtual TF1 * Get(size_t nFunc)
Returns a function sum of nFunc base functions.
ICARUShitFitCache fFitCache
Cached functions for multi-peak fits.
std::vector< ICARUSPeakFitParams_t > ICARUSPeakParamsVec
std::string to_string(WindowPattern const &pattern)
double ComputeChiSquare(TF1 func, TH1 *histo) const
void hit::ICARUSHitFinder::produce ( art::Event &  evt)

Definition at line 322 of file ICARUSHitFinder_module.cc.

323  {
324 
325  //0
326  //return;
327 
328  std::ofstream output("areaFit.out");
329 
330  // std::cout << " ICARUSHitFinder produce " << std::endl;
331 
332  //GET THE GEOMETRY.
333  art::ServiceHandle<geo::Geometry> geom;
334 
335  // ###############################################
336  // ### Making a ptr vector to put on the event ###
337  // ###############################################
338  // this contains the hit collection
339  // and its associations to wires and raw digits
340 
341  // Handle the filtered hits collection...
343 
344  // if (fAllHitsInstanceName != "") filteredHitCol = &hcol;
345 
346  // ##########################################
347  // ### Reading in the Wire List object(s) ###
348  // ##########################################
349  art::Handle< std::vector<recob::Wire> > wireVecHandle;
350  evt.getByLabel(fCalDataModuleLabel,wireVecHandle);
351 
352  // #################################################################
353  // ### Reading in the RawDigit associated with these wires, too ###
354  // #################################################################
355  art::FindOneP<raw::RawDigit> RawDigits
356  (wireVecHandle, evt, fCalDataModuleLabel);
357 
358  // Channel Number
360 
361 
362  std::vector<float> holder; //HOLDS SIGNAL DATA.
363  std::vector<short> rawadc; //UNCOMPRESSED ADC VALUES.
364 
365  std::vector<float> startTimes; //STORES TIME OF WINDOW START.
366  std::vector<float> maxTimes; //STORES TIME OF LOCAL MAXIMUM.
367  std::vector<float> endTimes; //STORES TIME OF WINDOW END.
368  std::vector<float> peakHeight; //STORES ADC COUNT AT THE MAXIMUM.
369  std::vector<float> hitrms; //STORES CHARGE WEIGHTED RMS OF TIME ACROSS THE HIT.
370  std::vector<double> charge; //STORES THE TOTAL CHARGE ASSOCIATED WITH THE HIT.
371 
372  // geo::SigType_t sigType = geo::kInduction;
373  std::stringstream numConv;
374 
375  unsigned int hwC(0),lwC(5728); //lowest and highest wire with physical deposition in Collection
376  unsigned int hwI2(0),lwI2(5728); //same for Induction
377  unsigned int hwI1(0),lwI1(2112);
378 
379  unsigned int nw1hitC(0),nw1hitI2(0),nw1hitI1(0); //number of wires (between lowest and highest) with a single hit (rough definition of hitfinding efficiency)
380  int nhWire[5600];
381  for(int jw=0;jw<5600;jw++)
382  nhWire[jw]=0;
383  float wCharge[5600];
384  for(int jw=0;jw<5600;jw++)
385  wCharge[jw]=0;
386  float wInt[5600];
387  for(int jw=0;jw<5600;jw++)
388  wInt[jw]=0;
389 
390 
391  unsigned int nhitsC(0),nhitsI1(0),nhitsI2(0); //total number of reconstructed hits in a view
392 
393  unsigned int nnhitsC(0),nnhitsI1(0),nnhitsI2(0); //total number of reconstructed hits in a view
394 
395 
396  unsigned int minWireC,maxWireC;
397  if(fThetaAngle==45) {minWireC=2539; maxWireC=3142;}
398  if(fThetaAngle==0) {minWireC=2535; maxWireC=4486;}
399  if(fThetaAngle==20) {minWireC=2539; maxWireC=4000;}
400  if(fThetaAngle==40) {minWireC=2539; maxWireC=3190;}
401  if(fThetaAngle==60) {minWireC=2539; maxWireC=2905;}
402  if(fThetaAngle==70) {minWireC=2539; maxWireC=2805;}
403  if(fThetaAngle==80) {minWireC=2539; maxWireC=2740;}
404 
405  //GET THE LIST OF BAD CHANNELS.
406  lariov::ChannelStatusProvider const& channelStatus
407  = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
408 
410  = channelStatus.BadChannels();
411 
412  unsigned int minWireI2=2539; //empirical
413  unsigned int maxWireI2=4700;
414  unsigned int minDrift=850;
415  unsigned int maxDrift=1500;
416 
417  //### Looping over the wires ###
418  //##############################
419  for(size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++)
420  {
421  // ####################################
422  // ### Getting this particular wire ###
423  // ####################################
424  art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
425  art::Ptr<raw::RawDigit> rawdigits = RawDigits.at(wireIter);
426 
427  // --- Setting Channel Number and Signal type ---
428  channel = wire->Channel();
429 
430  std::vector<float> signal(wire->Signal());
431 
432 
433  // get the WireID for this hit
434  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
435  // for now, just take the first option returned from ChannelToWire
436  geo::WireID wid = wids[0];
437  // We need to know the plane to look up parameters
438  geo::PlaneID::PlaneID_t plane = wid.Plane;
439  size_t cryostat=wid.Cryostat;
440  size_t tpc=wid.TPC;
441  size_t iWire=wid.Wire;
442 
443 
444  holder.clear();
445  localmeans.clear();
446 
447  //GET THE REFERENCE TO THE CURRENT raw::RawDigit.
448  channel = rawdigits->Channel();
449  fDataSize = rawdigits->Samples();
450 
451  std::vector<geo::WireID> widVec = geom->ChannelToWire(channel);
452  size_t iwire = widVec[0].Wire;
453  // size_t plane = widVec[0].Plane;
454 
455 
456  rawadc.resize(fDataSize);
457  holder.resize(fDataSize);
458 
459  //UNCOMPRESS THE DATA.
460  if (fUncompressWithPed) {
461  int pedestal = (int)rawdigits->GetPedestal();
462  raw::Uncompress(rawdigits->ADCs(), rawadc, pedestal, rawdigits->Compression());
463  }
464  else{
465  raw::Uncompress(rawdigits->ADCs(), rawadc, rawdigits->Compression());
466  }
467 
468  mf::LogDebug("ICARUSHitFinder") << " pedestal " <<rawdigits->GetPedestal() << std::endl;
469 
470  for(unsigned int bin = 0; bin < fDataSize; ++bin){
471  //holder[bin]=(rawadc[bin]-rawdigits->GetPedestal());
472  holder[bin]=signal[bin];
473 
474  //std::cout << " bin " << bin << " rawadc " << rawadc[bin]-rawdigits->GetPedestal() << " signal " << signal[bin] << std::endl;
475  if(plane == 0) holder[bin]=-holder[bin];
476  //if(plane == 1) holder[bin]=-holder[bin];
477  // if(plane==2&&iwire==2600) std::cout << " wire 2600 bin " << bin << " signal " << holder[bin] << std::endl;
478  }
479 
480  if(plane==0&&iwire<lwI1) lwI1=iwire;
481  if(plane==0&&iwire>hwI1) hwI1=iwire;
482  if(plane==1&&iwire<lwI2) lwI2=iwire;
483  if(plane==1&&iwire>hwI2) hwI2=iwire;
484  if(plane==2&&iwire<lwC) lwC=iwire;
485  if(plane==2&&iwire>hwC) hwC=iwire;
486 
487 
488 
489  // sigType = geom->SignalType(channel);
490 
491  peakHeight.clear();
492  endTimes.clear();
493  startTimes.clear();
494  maxTimes.clear();
495  charge.clear();
496  hitrms.clear();
497 
498  bool channelSwitch = false;
499 
500  for(auto it = BadChannels.begin(); it != BadChannels.end(); it++)
501  {
502  if(channel==*it)
503  {
504  channelSwitch = true;
505  break;
506  }
507  }
508 
509  if(channelSwitch==false)
510  { //1
511 
512  // Hit finding parameters
513  double chargeErr(0); //CHI2/NDF and error on charge.
514  // double hrms(0);
515  //double totSig;
516 
517  double chi2null=ComputeNullChiSquare(holder);
518 // std::cout << " wire " << iWire << " chi2null " << chi2null << std::endl;
519  fNullChi2->Fill(chi2null);
520 
521 
522 
524 
526 
527  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
528 
529  std::vector<float> tempVec = holder;
530  recob::Wire::RegionsOfInterest_t::datarange_t rangeData(size_t(0),std::move(tempVec));
531 
532  fHitFinderTool->findHitCandidates(rangeData, 0,channel,0,hitCandidateVec);
533  //int jc=0;
534  for(auto& hitCand : hitCandidateVec) {
535  expandHit(hitCand,holder,hitCandidateVec);
536 
537  }
538 
539 
540 
541  fHitFinderTool->MergeHitCandidates(rangeData, hitCandidateVec, mergedCandidateHitVec);
542 
543  //numHits = hits.size();
544  int nghC=0;
545  int nghI2=0;
546  int nghI1=0;
547 
548  // if(plane==0)
549  // std::cout << " plane " << plane << " Wire " << iwire << " numhits " << hitCandidateVec.size() <<" mergedhits " << mergedCandidateHitVec.size() << std::endl;
550 
551 
552  //FIT ONLY COLLECTION HITS
553  if(plane==2) {
554  //std::cout << " mergedcands size " << mergedCandidateHitVec.size() << std::endl;
555 
556  for(auto& mergedCands : mergedCandidateHitVec)
557  {
558 
559 
560  int startT= mergedCands.front().startTick-fFittingRange;
561  int endT = mergedCands.back().stopTick+fFittingRange;
562  //std::cout << " fitting range " << fFittingRange << std::endl;
563 
564  float mean;
565  computeBestLocalMean(mergedCands,holder,mergedCandidateHitVec,mean);
566  localmeans.push_back(mean);
567  mf::LogDebug("ICARUSHitFinder") << " adding localmean " << mean << std::endl;
568  // ### Putting in a protection in case things went wrong ###
569  // ### In the end, this primarily catches the case where ###
570  // ### a fake pulse is at the start of the ROI ###
571  if (endT - startT < 5) continue;
572  fWidthC->Fill(endT-startT);
573  // #######################################################
574  // ### Clearing the parameter vector for the new pulse ###
575  // #######################################################
576 
577  // === Setting the number of Gaussians to try ===
578  int nGausForFit = mergedCands.size();
579 
580  // ##################################################
581  // ### Calling the function for fitting ICARUS ###
582  // ##################################################
583  double chi2PerNDF(0.);
584  double chi2Long(0.);
585 
586  int NDF(1);
587  const int npk=mergedCands.size();
588  ICARUSPeakParamsVec peakParamsVec(npk);
589  peakParamsVec.clear();
590  int islong=0;
591  if (mergedCands.size() <= fMaxMultiHit)
592  {
593 
594  // std::cout << "setting iwire " << iwire << std::endl;
595  // fPeakFitterTool->setWire(iwire);
596  // std::cout << " fitting iwire " << iwire << std::endl;
597  // std::cout << " cryostat " << cryostat << " tpc " << tpc << " plane " << plane << " wire " << iwire << std::endl;
598  findMultiPeakParameters(signal, mergedCands, peakParamsVec, chi2PerNDF, NDF, iwire);
599 
600  if (!(chi2PerNDF < std::numeric_limits<double>::infinity()))
601  {
602  chi2PerNDF = 200.;
603  NDF = 2;
604  }
605  fFirstChi2->Fill(chi2PerNDF);
606  // std::cout << " wire " << iwire << " first chi2NDF " << chi2PerNDF << std::endl;
607  }
608 
609  // std::cout << " before longpulse chi2 " << chi2PerNDF << " threshold " << fChi2NDF << std::endl;
610  if (chi2PerNDF < 10.)
611  fChi2->Fill(chi2PerNDF);
612  // if(chi2PerNDF<0.1) // change from 10 to reduce output
613  // std::cout << " wire " << iwire << " SMALL chi2NDF " << chi2PerNDF << " thr chi2NDF " << fChi2NDF << std::endl;
614  ICARUSPeakParamsVec peakParamsLong(npk);
615  peakParamsLong.clear();
616  if (chi2PerNDF > fChi2NDF)
617  {
618  islong=1;
619  findLongPeakParameters(signal, mergedCands, peakParamsLong, chi2Long, NDF, iwire);
620  // if(chi2Long<0.3) std::cout << " small chi2long " << chi2Long << std::endl;
621  if(chi2Long<chi2PerNDF&&chi2Long>0.1) {
622  fChi2->Fill(chi2Long);
623  peakParamsVec=peakParamsLong;
624  }
625  else { fChi2->Fill(chi2PerNDF);
626  // if(chi2PerNDF<0.1) // change from 10 to reduce output
627  // std::cout << " wire " << iwire << " SMALLSMALL chi2NDF " << chi2PerNDF << " thr chi2NDF " << fChi2NDF << std::endl;
628  }
629  }
630 
631  // unsigned int jhit=0;
632 
633  //std::cout << " before peak loop" << std::endl;
634  // for(const auto& peakParams : peakParamsVec)
635 for(unsigned int jhit=0;jhit<mergedCands.size(); jhit++)
636  {
637  //float fitCharge=chargeFunc(peakMean, peakAmp, peakWidth, fAreaNormsVec[plane],startT,endT);
638  //float fitChargeErr = std::sqrt(TMath::Pi()) * (peakAmpErr*peakWidthErr + peakWidthErr*peakAmpErr);
639  unsigned int startInt=mergedCands[jhit].startTick-fIntegratingRange;
640  unsigned int endInt=mergedCands[jhit].stopTick+fIntegratingRange;
641 
642  if(jhit>=1&&startInt<mergedCands[jhit-1].stopTick) startInt=mergedCands[jhit-1].stopTick;
643  if(jhit<mergedCands.size()-1&&endInt>mergedCands[jhit+1].startTick) endInt=mergedCands[jhit+1].startTick;
644 
645  float fitCharge=0;
646  float peakAmp, peakMean, peakLeft, peakRight, peakBaseline;
647  float peakSlope=0, peakFitWidth=0;
648  float peakMeanErr, peakAmpErr;
649  if(!islong) {
650  // TF1 Func("ICARUSfunc",fitf,start,end,1+5*mergedCands.size());
651  TF1& Func = *(fFitCache.Get(mergedCands.size()));
652  assert(&Func);
653  Func.SetParameter(0, mergedCands.size());
654 
655  float intBaseline=0;
656 
657  ICARUSPeakFitParams_t peakParams=peakParamsVec[jhit];
658 
659  // Extract values for this FITTED peak
660  peakAmp = peakParams.peakAmplitude;
661  peakMean = peakParams.peakCenter;
662  //float peakWidth = peakParams.peakSigma;
663  peakLeft = peakParams.peakTauLeft;
664  peakRight = peakParams.peakTauRight;
665  peakBaseline = peakParams.peakBaseline;
666 
667 
668  // Place one bit of protection here
669  if (std::isnan(peakAmp))
670  {
671  // std::cout << "**** hit peak amplitude is a nan! Channel: " << channel << ", start tick: " << startT << std::endl;
672  continue;
673  }
674 
675  // Extract errors
676  peakAmpErr = peakParams.peakAmplitudeError;
677  peakMeanErr = peakParams.peakCenterError;
678  // float peakWidthErr = peakParams.peakSigmaError;
679  Func.SetParameter(1+5*jhit,peakBaseline);
680  Func.SetParameter(2+5*jhit,peakAmp);
681  Func.SetParameter(3+5*jhit,peakMean);
682  Func.SetParameter(4+5*jhit,peakRight);
683  Func.SetParameter(5+5*jhit,peakLeft);
684 
685  intBaseline+=(endInt-startInt)*peakBaseline;
686 
687  try
688  {
689  fitCharge=Func.Integral(startInt,endInt)-(endInt-startInt)*localmeans[jhit];
690 
691  }
692  catch(...) {
693  mf::LogWarning("ICARUSHitFinder") << "Icarus numerical 32 failed";
694  fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.)-(endInt-startInt)*localmeans[jhit];
695  }
696  }
697  else {
698  // TF1 FuncLong("ICARUSfuncLong",fitlong,start,end,1+7*mergedCands.size());
699  TF1& FuncLong = *(fLongFitCache.Get(mergedCands.size()));
700  assert(&FuncLong);
701  FuncLong.SetParameter(0, mergedCands.size());
702 
703 
704  float intBaseline=0;
705 
706  ICARUSPeakFitParams_t peakParams=peakParamsVec[jhit];
707 
708  // Extract values for this FITTED peak
709  peakAmp = peakParams.peakAmplitude;
710  peakMean = peakParams.peakCenter;
711  //float peakWidth = peakParams.peakSigma;
712  peakLeft = peakParams.peakTauLeft;
713  peakRight = peakParams.peakTauRight;
714  peakBaseline = peakParams.peakBaseline;
715  intBaseline+=(endInt-startInt)*peakBaseline;
716  peakAmpErr = peakParams.peakAmplitudeError;
717  peakMeanErr = peakParams.peakCenterError;
718 
719  FuncLong.SetParameter(1+7*jhit,peakBaseline);
720  FuncLong.SetParameter(2+7*jhit,peakAmp);
721  FuncLong.SetParameter(3+7*jhit,peakMean);
722  FuncLong.SetParameter(4+7*jhit,peakRight);
723  FuncLong.SetParameter(5+7*jhit,peakLeft);
724  FuncLong.SetParameter(6+7*jhit,peakFitWidth);
725  FuncLong.SetParameter(7+7*jhit,peakSlope);
726 
727 
728  try
729  { fitCharge=FuncLong.Integral(startInt,endInt)-(endInt-startInt)*localmeans[jhit];
730  //fitCharge=FuncLong.Integral(start-35,end+35);
731  }
732  catch(...)
733  {mf::LogWarning("ICARUSHitFinder") << "Icarus numerical integration failed";
734  fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.)-(endInt-startInt)*localmeans[jhit];
735  }
736  }
737  if(isnan(fitCharge)&&!islong) fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.);
738  if(isnan(fitCharge)&&islong) fitCharge=std::accumulate(holder.begin() + (int) startInt, holder.begin() + (int) endInt, 0.);
739  //Func.Integral(start,end);
740 
741  // float totSig20=std::accumulate(holder.begin() + (int) start-35, holder.begin() + (int) end+35, 0.);
742 
743  float totSig=std::accumulate(holder.begin()+ (int) startInt, holder.begin()+ (int) endInt, 0.)-(endInt-startInt)*localmeans[jhit];
744  fBaselineC->Fill(localmeans[jhit]);
745  // float fitChargeErr=0;
746  // if(plane==2&&iWire==2842)
747 
748 
749 
750 //std::cout << " before hit creator " << std::endl;
752  *wire, //RAW DIGIT REFERENCE.
753  wid, //WIRE ID.
754  startInt, //START TICK.
755  endInt, //END TICK.
756  (peakLeft+peakRight)/2., //RMS.
757  peakMean, //PEAK_TIME.
758  peakMeanErr, //SIGMA_PEAK_TIME.
759  peakAmp, //PEAK_AMPLITUDE.
760  peakAmpErr, //SIGMA_PEAK_AMPLITUDE.
761  fitCharge, //HIT_INTEGRAL.
762  chargeErr, //HIT_SIGMA_INTEGRAL.
763  totSig, //SUMMED CHARGE.
764  nGausForFit, //MULTIPLICITY.
765  jhit, //LOCAL_INDEX.
766  chi2PerNDF, //WIRE ID.
767  NDF //DEGREES OF FREEDOM.
768  );
769 
770  mf::LogDebug("ICARUSHitFinder") << " fitcharge " << fitCharge << " totSig " << totSig << std::endl;
771  //std::cout << " summedADC " << hit.summedADC() << " integral " << hit.Integral() << std::endl;
772 
773  // filteredHitVec.push_back(hit.copy());
774  hcol.emplace_back(hit.move(), wire, rawdigits);
775 
776  bool wireWindowC=iWire>=minWireC&&iWire<=maxWireC;
777  bool outWireWindowC=iWire<=minWireC||iWire>=maxWireC;
778 
779  if(wireWindowC) {
780  nghC++;
781  if(nghC==1) nw1hitC++;
782  nhWire[iWire]++;
783  fHeightC->Fill(peakAmp);
784  // fWidthC->Fill(peakWidth);
785  // fAreaC->Fill(totSig);
786  wCharge[iWire]+=fitCharge;
787  wInt[iWire]+=totSig;
788  }
789 
790 
791  nhitsC++;
792 
793 
794  if(tpc==0&&cryostat==0&&outWireWindowC) {
795  nnhitsC++; fNoiseC->Fill(peakAmp); }
796 
797 
798 
799  if(cryostat==0&&tpc==0)
800  if(wCharge[iwire]>0)
801  fAreaC->Fill(wCharge[iwire]);
802  if(cryostat==0&&tpc==0)
803  if(wInt[iwire]>0)
804  fIntegralC->Fill(wInt[iwire]);
805  if(cryostat==0&&tpc==0)
806  if(wCharge[iwire]>0)
807  output << iwire << " " <<wInt[iwire] << std::endl;
808  if(wInt[iwire]>0&&cryostat==0&&tpc==0)
809  fAreaInt->Fill(wCharge[iwire]/wInt[iwire]);
810 
811  } // loop on peakparams vector
812 
813  } // loop on merged hits
814 
815  // std::cout << " after coll hit " << std::endl;
816  } //COLLECTION
817 
818  if(plane==0||plane==1) {
819  for(auto& mergedCands : mergedCandidateHitVec)
820  {
821  for(size_t jh=0;jh<mergedCands.size();jh++) {
822  // std::cout << " plane " << plane << " wire " << iwire << " hit " <<mergedCands[jh].hitCenter << std::endl;
823  //FOR INDUCTION HITS STORE RAW INFORMATION
825  *wire, //RAW DIGIT REFERENCE.
826  wid, //WIRE ID.
827  mergedCands[jh].startTick, //START TICK.
828  mergedCands[jh].stopTick, //END TICK.
829  mergedCands[jh].hitSigma, //RMS.
830  mergedCands[jh].hitCenter, //PEAK_TIME.
831  0, //SIGMA_PEAK_TIME.
832  mergedCands[jh].hitHeight, //PEAK_AMPLITUDE.
833  0, //SIGMA_PEAK_AMPLITUDE.
834  std::accumulate(holder.begin() + (int) mergedCands[jh].startTick, holder.begin() + (int) mergedCands[jh].stopTick, 0.), //HIT_INTEGRAL.
835  0, //HIT_SIGMA_INTEGRAL.
836  std::accumulate(holder.begin() + (int) mergedCands[jh].startTick, holder.begin() + (int) mergedCands[jh].stopTick, 0.), //SUMMED CHARGE.
837  1, //MULTIPLICITY.
838  jh, //LOCAL_INDEX.
839  0, //WIRE ID.
840  int(mergedCands[jh].stopTick-mergedCands[jh].startTick+1) //DEGREES OF FREEDOM.
841  );
842  hcol.emplace_back(hit.move(), wire, rawdigits);
843 
844  bool driftWindow=(mergedCands[jh].hitCenter)>=minDrift&&(mergedCands[jh].hitCenter)<=maxDrift;
845  bool wireWindowI2=iWire>=minWireI2&&iWire<=maxWireI2;
846  bool outDriftWindow=(mergedCands[jh].hitCenter)<=minDrift||(mergedCands[jh].hitCenter)>=maxDrift;
847  bool outWireWindowI2=iWire<=minWireI2||iWire>=maxWireI2;
848 
849 
850  if(plane==0&&driftWindow) {
851  // std::cout << " wire " << hits[i].iWire << " ngh " << ngh << std::endl;
852  nghI1++;
853  if(nghI1==1) nw1hitI1++;
854 
855  // std::cout << " wire " << hits[i].iWire << " nw1h " << nw1 << std::endl;
856  fHeightI1->Fill(mergedCands[jh].hitHeight);
857  fWidthI1->Fill(mergedCands[jh].hitSigma);
858  }
859  if(plane==1&&wireWindowI2) {
860  // std::cout << " wire " << hits[i].iWire << " ngh " << ngh << std::endl;
861  nghI2++;
862  if(nghI2==1) nw1hitI2++;
863  // std::cout << " filling height histo wire" << hits[i].iWire << " drift " << hits[i].hitCenter << " amplitude " << amplitude << std::endl;
864  fHeightI2->Fill(mergedCands[jh].hitHeight);
865  fWidthI2->Fill(mergedCands[jh].hitSigma);
866  }
867 
868 
869  if(plane==0) nhitsI1++;
870  if(plane==1) nhitsI2++;
871 
872  if(plane==0&&tpc==0&&cryostat==0&&outDriftWindow) {nnhitsI1++; fNoiseI1->Fill(mergedCands[jh].hitHeight); }
873  if(plane==1&&tpc==0&&cryostat==0&&outWireWindowI2) { nnhitsI2++; fNoiseI2->Fill(mergedCands[jh].hitHeight); }
874  }} // merged loop
875  } //INDUCTION
876 
877 
878 
879  } //end channel condition
880 
881  } //end loop on channels
882 
883 // std::cout << " nhitsI1 " << nhitsI1 <<" nhitsI2 " << nhitsI2 <<" nhitsC " << nhitsC << std::endl;
884 
885 
886  for(unsigned int jw=minWireC;jw<maxWireC;jw++)
887  fnhwC->Fill(nhWire[jw]);
888 
889 
890  hcol.put_into(evt);
891  //std::cout << " end ICARUSHitfinder " << std::endl;
892 
893 
894  } //end produce
size_t fMaxMultiHit
maximum hits for multi fit
std::set< raw::ChannelID_t > ChannelSet_t
Type of set of channel IDs.
unsigned int PlaneID_t
Type for the ID number.
Definition: geo_types.h:473
struct ICARUSPeakFitParams{float peakCenter ICARUSPeakFitParams_t
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
process_name hit
Definition: cheaterreco.fcl:51
virtual ChannelSet_t BadChannels() const =0
Returns a copy of set of bad channel IDs for the current run.
void computeBestLocalMean(std::vector< reco_tool::ICandidateHitFinder::HitCandidate > h, std::vector< float > holder, reco_tool::ICandidateHitFinder::MergeHitCandidateVec how, float &localmean)
constexpr details::BinObj< T > bin(T value)
Returns a wrapper to print the specified data in binary format.
double ComputeNullChiSquare(std::vector< float >) const
ICARUSlongHitFitCache fLongFitCache
Cached functions for long hits.
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:83
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:32
A class handling a collection of hits and its associations.
Definition: HitCreator.h:508
int fFittingRange
semi-width of interval where to fit hit
std::unique_ptr< reco_tool::ICandidateHitFinder > fHitFinderTool
For finding candidate hits.
std::vector< double > localmeans
Class providing information about the quality of channels.
virtual TF1 * Get(size_t nFunc)
Returns a function sum of nFunc base functions.
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
ICARUShitFitCache fFitCache
Cached functions for multi-peak fits.
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:13
std::vector< ICARUSPeakFitParams_t > ICARUSPeakParamsVec
BEGIN_PROLOG sequence::SlidingWindowTriggerPatternsOppositeWindows END_PROLOG simSlidingORM6O6 effSlidingORW output
void expandHit(reco_tool::ICandidateHitFinder::HitCandidate &h, std::vector< float > holder, std::vector< reco_tool::ICandidateHitFinder::HitCandidate > how)
TCEvent evt
Definition: DataStructs.cxx:8
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
void findLongPeakParameters(const std::vector< float > &, const reco_tool::ICandidateHitFinder::HitCandidateVec &, ICARUSPeakParamsVec &, double &, int &, int) const
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
std::vector< HitCandidateVec > MergeHitCandidateVec
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
Definition: raw.cxx:776
int fIntegratingRange
semi-width of interval where to integrate fitting function
void findMultiPeakParameters(const std::vector< float > &, const reco_tool::ICandidateHitFinder::HitCandidateVec &, ICARUSPeakParamsVec &, double &, int &, int) const
std::vector< HitCandidate > HitCandidateVec
void hit::ICARUSHitFinder::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 249 of file ICARUSHitFinder_module.cc.

250  {
251  fUncompressWithPed = p.get< bool >("UncompressWithPed", true);
252  fDigitModuleLabel = p.get< art::InputTag >("DigitModuleLabel", "daq");
253  fCalDataModuleLabel = p.get< std::string >("CalDataModuleLabel");
254  fMaxMultiHit = p.get< int >("MaxMultiHit");
255  fChi2NDF = p.get< double >("Chi2NDF");
256  fLongPulseWidthVec = p.get< std::vector<int>>("LongPulseWidth", std::vector<int>() = {16,16,16});
257  fLongMaxHitsVec = p.get< std::vector<int>>("LongMaxHits", std::vector<int>() = {25,25,25});
258  fThetaAngle=p.get< double >("ThetaAngle");
259  fMinWidth=p.get< double >("MinWidth");
260  fMaxWidthMult=p.get< double >("MaxWidthMult");
261  fFittingRange=p.get< int >("FittingRange");
262  fIntegratingRange=p.get< int >("IntegratingRange");
263 
264 
265  fHitFinderTool = art::make_tool<reco_tool::ICandidateHitFinder>(p.get<fhicl::ParameterSet>("CandidateHits"));
266  // fPeakFitterTool = art::make_tool<reco_tool::PeakFitterICARUS>(p.get<fhicl::ParameterSet>("PeakFitter"));
267 
268  mf::LogInfo("ICARUSHitFinder_module") << "fDigitModuleLabel: " << fDigitModuleLabel << std::endl;
269  }
size_t fMaxMultiHit
maximum hits for multi fit
std::vector< int > fLongPulseWidthVec
Sets width of hits used to describe long pulses.
pdgs p
Definition: selectors.fcl:22
double fMinWidth
minimum initial width for ICARUS fit
double fMaxWidthMult
multiplier for max width for ICARUS fit
int fFittingRange
semi-width of interval where to fit hit
std::unique_ptr< reco_tool::ICandidateHitFinder > fHitFinderTool
For finding candidate hits.
std::vector< int > fLongMaxHitsVec
maximum Chisquared / NDF allowed for a hit to be saved
int fIntegratingRange
semi-width of interval where to integrate fitting function
void hit::ICARUSHitFinder::setWire ( int  i)
inline

Definition at line 171 of file ICARUSHitFinder_module.cc.

171  {
172  iWire=i;
173  // std::cout << " setting iwire " << iWire << std::endl;
174  } ;

Member Data Documentation

TH1F* hit::ICARUSHitFinder::fAreaC
private

Definition at line 204 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fAreaInt
private

Definition at line 207 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fBaselineC
private

Definition at line 208 of file ICARUSHitFinder_module.cc.

std::string hit::ICARUSHitFinder::fCalDataModuleLabel
private

Definition at line 182 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fChi2
private

Definition at line 200 of file ICARUSHitFinder_module.cc.

double hit::ICARUSHitFinder::fChi2NDF
private

Definition at line 192 of file ICARUSHitFinder_module.cc.

unsigned int hit::ICARUSHitFinder::fDataSize
private

Definition at line 177 of file ICARUSHitFinder_module.cc.

art::InputTag hit::ICARUSHitFinder::fDigitModuleLabel
private

Definition at line 178 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fFirstChi2
private

Definition at line 197 of file ICARUSHitFinder_module.cc.

ICARUShitFitCache hit::ICARUSHitFinder::fFitCache
mutableprivate

Cached functions for multi-peak fits.

Definition at line 228 of file ICARUSHitFinder_module.cc.

int hit::ICARUSHitFinder::fFittingRange
private

semi-width of interval where to fit hit

Definition at line 222 of file ICARUSHitFinder_module.cc.

const geo::GeometryCore* hit::ICARUSHitFinder::fGeometry = lar::providerFrom<geo::Geometry>()
private

Definition at line 231 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fHeightC
private

Definition at line 201 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fHeightI1
private

Definition at line 215 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fHeightI2
private

Definition at line 212 of file ICARUSHitFinder_module.cc.

std::unique_ptr<reco_tool::ICandidateHitFinder> hit::ICARUSHitFinder::fHitFinderTool
private

For finding candidate hits.

Definition at line 188 of file ICARUSHitFinder_module.cc.

std::string hit::ICARUSHitFinder::fHitLabelName
private

Definition at line 183 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fIntegralC
private

Definition at line 206 of file ICARUSHitFinder_module.cc.

int hit::ICARUSHitFinder::fIntegratingRange
private

semi-width of interval where to integrate fitting function

Definition at line 223 of file ICARUSHitFinder_module.cc.

ICARUSlongHitFitCache hit::ICARUSHitFinder::fLongFitCache
mutableprivate

Cached functions for long hits.

Definition at line 229 of file ICARUSHitFinder_module.cc.

std::vector<int> hit::ICARUSHitFinder::fLongMaxHitsVec
private

maximum Chisquared / NDF allowed for a hit to be saved

Maximum number hits on a really long pulse train

Definition at line 193 of file ICARUSHitFinder_module.cc.

std::vector<int> hit::ICARUSHitFinder::fLongPulseWidthVec
private

Sets width of hits used to describe long pulses.

Definition at line 194 of file ICARUSHitFinder_module.cc.

size_t hit::ICARUSHitFinder::fMaxMultiHit
private

maximum hits for multi fit

Definition at line 191 of file ICARUSHitFinder_module.cc.

double hit::ICARUSHitFinder::fMaxWidthMult
private

multiplier for max width for ICARUS fit

Definition at line 221 of file ICARUSHitFinder_module.cc.

double hit::ICARUSHitFinder::fMinWidth
private

minimum initial width for ICARUS fit

Definition at line 220 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fnhwC
private

Definition at line 205 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fNoiseC
private

Definition at line 203 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fNoiseI1
private

Definition at line 217 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fNoiseI2
private

Definition at line 214 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fNullChi2
private

Definition at line 198 of file ICARUSHitFinder_module.cc.

std::string hit::ICARUSHitFinder::fSpillName
private

Definition at line 179 of file ICARUSHitFinder_module.cc.

int hit::ICARUSHitFinder::fThetaAngle
private

Definition at line 185 of file ICARUSHitFinder_module.cc.

bool hit::ICARUSHitFinder::fUncompressWithPed
private

Definition at line 186 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fWidthC
private

Definition at line 202 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fWidthI1
private

Definition at line 216 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fWidthI2
private

Definition at line 213 of file ICARUSHitFinder_module.cc.

TH1F* hit::ICARUSHitFinder::fWire2771
private

Definition at line 210 of file ICARUSHitFinder_module.cc.

int hit::ICARUSHitFinder::iWire
private

Definition at line 226 of file ICARUSHitFinder_module.cc.

std::vector<double> hit::ICARUSHitFinder::localmeans
private

Definition at line 174 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakAmplitude

Definition at line 143 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakAmplitudeError

Definition at line 144 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakBaseline

Definition at line 153 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakBaselineError

Definition at line 154 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakCenterError

Definition at line 140 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakFitWidth

Definition at line 149 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakFitWidthError

Definition at line 150 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakSigma

Definition at line 141 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakSigmaError

Definition at line 142 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakSlope

Definition at line 151 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakSlopeError

Definition at line 152 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakTauLeft

Definition at line 145 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakTauLeftError

Definition at line 146 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakTauRight

Definition at line 147 of file ICARUSHitFinder_module.cc.

float hit::ICARUSHitFinder::peakTauRightError

Definition at line 148 of file ICARUSHitFinder_module.cc.


The documentation for this class was generated from the following file: