TimeBasedDisambig.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // TimeBasedDisambig.cxx
4 //
5 // trj@fnal.gov
6 // tjyang@fnal.gov
7 // iamjaejunkim@gmail.com
8 //
9 // description
10 //
11 // Based on time based hit matching algorithm
12 //
13 ////////////////////////////////////////////////////////////////////////
14 
15 
16 //Framework includes:
19 #include "art_root_io/TFileService.h"
20 #include "art_root_io/TFileDirectory.h"
22 
30 #include "TimeBasedDisambig.h"
31 //#include "MCCheater/BackTracker.h"
33 
34 #include <map>
35 #include <cmath>
36 #include <vector>
37 #include <iostream>
38 #include <fstream>
39 #include <cstdlib>
40 
41 #include "TH1D.h"
42 #include "TF1.h"
43 
44 using namespace std;
45 
46 namespace dune{
47 
48 TimeBasedDisambig::TimeBasedDisambig(fhicl::ParameterSet const& pset)
49 {
50  this->reconfigure(pset);
51 }
52 
53 //----------------------------------------------------------
55 {
56 
57  fTimeCut = p.get<double>("TimeCut");
58  fDistanceCut = p.get<double>("DistanceCut");
59  fDistanceCutClu = p.get<double>("DistanceCutClu");
60 }
61 
62 
63 //----------------------------------------------------------
64 //----------------------------------------------------------
65 void TimeBasedDisambig::RunDisambig( const std::vector< art::Ptr<recob::Hit> > &OrigHits )
66 //void TimeBasedDisambig::RunDisambig()
67 {
68  //fDisambigHits.clear();
69 
70  //create geometry and backtracker servicehandle object
72  //art::ServiceHandle<cheat::BackTracker> bt;
73  //define timeoffset and the maximum chargeratio allowed when matching induction plane hits
74  //the timeoffset array is going to be replaced with defined values in the detectorproperties package later
75  double timeoffset[3]={10.4,5.2,0.0};
76  double rmax=2.5;
77 
78  //create vectors for induction plane and collection plane hits
79  std::vector<art::Ptr<recob::Hit> > hitsUV;
80  std::vector<art::Ptr<recob::Hit> > hitsZ;
81 
82  //create a vector where the disambiguated hits are going to be stored
83  //std::vector<art::Ptr<recob::Hit> > DisambiguatedHits;
84  //std::unique_ptr<std::vector<recob::Hit> > DisambiguatedHits(new std::vector<recob::Hit>);
85  //std::vector< <recob::Hit> > DisambiguatedHits;
86  std::vector< recob::Hit > DisambiguatedHits;
87 
88  //double HitsTotal=0;
89  //double HitsCorrect=0;
90 
91  //loop over detector volumn in each component
92  for (unsigned int Cstat=0; Cstat < geo->Ncryostats(); ++Cstat){
93  for (unsigned int APA=0; APA < geo->Cryostat(Cstat).NTPC()/2; ++APA){
94  hitsUV.clear();
95  hitsZ.clear();
96  //loop over all induction and collection plane hits and save the hits within each component
97  for (size_t i = 0; i<OrigHits.size(); ++i){
98  if (OrigHits[i]->WireID().Cryostat==Cstat && (OrigHits[i]->WireID().TPC==APA*2 || OrigHits[i]->WireID().TPC==APA*2+1)){
99 
100  //save induction and collection plane hits in two separate vectors
101  if (OrigHits[i]->View()!=geo::kZ){
102  hitsUV.push_back(OrigHits[i]);
103  } else {
104  hitsZ.push_back(OrigHits[i]);
105  }
106 
107  }
108  }
109 
110  //define variables to save boundary information
111  double origin[3]={0.0,0.0,0.0};
112  double world[3]={0.0,0.0,0.0};
113  //double xbound[2]={-2000,2000};
114  //define xbound using peaktime for now but need to be replaced with coverted position in the future
115  double xbound[2]={0,2000};
116  double ybound[2]= {0.0,0.0};
117  double zbound[2]= {0.0,0.0};
118  geo->Cryostat(Cstat).TPC(APA*2).LocalToWorld(origin, world);
119  ybound[0]=world[1]-geo->Cryostat(Cstat).TPC(APA*2).HalfHeight();
120  ybound[1]=world[1]+geo->Cryostat(Cstat).TPC(APA*2).HalfHeight();
121  zbound[0]=world[2]-0.5*geo->Cryostat(Cstat).TPC(APA*2).Length();
122  zbound[1]=world[2]+0.5*geo->Cryostat(Cstat).TPC(APA*2).Length();
123 
124  //position-correction function needs to be applied on hits in certain time and longitudinal position ranges
125  //thus need to create a struct of vectors where hit informations including their positions returned in the
126  //disambiguation algorithm is going to be saved.
127  HitPos h;
128  std::vector<HitPos> InductionAll;
129  InductionAll.clear();
130 
131 
132  //loop over all induction plane hits
133  for (unsigned int uv0=0; uv0 < hitsUV.size(); ++uv0){
134  double PeakTimeMinusUV=hitsUV[uv0]->PeakTimePlusRMS(-1.)+timeoffset[1];
135  double PeakTimeUV=hitsUV[uv0]->PeakTime()+timeoffset[1];
136  if (hitsUV[uv0]->View()==geo::kU){
137  PeakTimeMinusUV+=timeoffset[1];
138  PeakTimeUV+=timeoffset[1];
139  }
140 
141  //define a varible to save charge per num
142  double ChargePerNum0=hitsUV[uv0]->Integral()/hitsUV[uv0]->Multiplicity();
143 
144  //loop over induction plane hits on a different view and among the induction plane hits whose charge per
145  //num to that of the hit in loop is within 2.5 find an induction plane hit whose peaktime is closest to
146  //that is in loop.
147  double MeanPeakTimeUV=99999;
148  unsigned int ChannelUV1=0;
149  for (unsigned int uv1=0; uv1 < hitsUV.size(); ++uv1){
150  double PeakTimeMinusUV1=hitsUV[uv1]->PeakTimePlusRMS(-1.)+timeoffset[1];
151  double PeakTimeUV1=hitsUV[uv1]->PeakTime()+timeoffset[1];
152  if (hitsUV[uv1]->View()==geo::kU){
153  PeakTimeMinusUV1+=timeoffset[1];
154  PeakTimeUV1+=timeoffset[1];
155  }
156  double rChargeUV=hitsUV[uv0]->Integral()/hitsUV[uv1]->Integral();
157  double ChargePerNum1=hitsUV[uv1]->Integral()/hitsUV[uv1]->Multiplicity();
158  if (hitsUV[uv0]->Multiplicity()>0 && hitsUV[uv1]->Integral()>0) rChargeUV=ChargePerNum0/ChargePerNum1;
159 
160  if (hitsUV[uv0]->View()!=hitsUV[uv1]->View() &&
161  fabs(PeakTimeMinusUV-PeakTimeMinusUV1)<MeanPeakTimeUV &&
162  fabs(PeakTimeUV-PeakTimeUV1)<20 &&
163  (rChargeUV<rmax && rChargeUV>(1/rmax)) ){
164  MeanPeakTimeUV=fabs(PeakTimeMinusUV-PeakTimeMinusUV1);
165  ChannelUV1=hitsUV[uv1]->Channel();
166  }
167  }
168 
169  //define variables to save positions returned
170  unsigned int ChannelUV0=hitsUV[uv0]->Channel();
171  double CandPosition[3]={0.0,0.0,0.0}; //candidate positions going to be returned in intersectionpoint function
172  double FinalPosition[3]={0.0,0.0,0.0}; //final position among the candidate positions
173  unsigned int FinalWireID=0;//final wireid for the induction plane hit in the loop
174  double MinVertBetUVandZ=99999;
175  //unsigned int FinalWire=0;
176 
177  //read wireid vectors of induction plane hits matched in the above loop
178  std::vector<geo::WireID> wireiduv0 = geo->ChannelToWire(ChannelUV0);
179  std::vector<geo::WireID> wireiduv1 = geo->ChannelToWire(ChannelUV1);
180 
181  //loop over the induction plane hit wireids, find candidate intersection positions using intersectionpoint
182  //function and find a position whose longitudinal position is closest to the average position returned above
183  //and save the positions and wireid in variables defined above
184  for (unsigned int wireid0=0; wireid0 < wireiduv0.size(); ++wireid0){
185  double VertPos=0; //vertical position of collection plane wire
186  double MinPeakTime=99999; //variable to store peaktime information temporarily
187  //loop over collection plane hits and among the collection plane hits whose charge per num to that
188  //of the hit in loop is within 2.5 find a collection plane hit that has closest time to the induction plane
189  //hit in the loop.
190  for (unsigned int z=0; z < hitsZ.size(); ++z){
191  double PeakTimeMinusZ=hitsZ[z]->PeakTimePlusRMS(-1.);
192  double PeakTimeZ=hitsZ[z]->PeakTime();
193  double ChargePerNumZ=hitsZ[z]->Integral()/hitsZ[z]->Multiplicity();
194  double rChargeZ=hitsUV[uv0]->Integral()/hitsZ[z]->Integral();
195  if (hitsUV[uv0]->Multiplicity()>0 && hitsZ[z]->Integral()>0) rChargeZ=ChargePerNum0/ChargePerNumZ;
196  double vert[3]={0.0,0.0,0.0};
197  std::vector<geo::WireID> wires = geo->ChannelToWire(hitsZ[z]->Channel());
198  geo->Cryostat(wires[0].Cryostat).TPC(wires[0].TPC).Plane(2).Wire(wires[0].Wire).GetCenter(vert);
199  if (fabs(PeakTimeMinusUV-PeakTimeMinusZ)<MinPeakTime && fabs(PeakTimeUV-PeakTimeZ)<20 &&
200  (rChargeZ<rmax && rChargeZ>(1/rmax)) &&
201  wireiduv0[wireid0].TPC==wires[0].TPC){
202  MinPeakTime=fabs(PeakTimeMinusUV-PeakTimeMinusZ);
203  VertPos=vert[2];
204  }
205  }
206 
207  //loop over collection plane hits and among the collection plane hits whose charge per num to that
208  //of the hit in the loop is within 2.5 find collection plane hits whose time difference to the induction
209  //plane hit is within 10 and take the average position of the collection plane hits.
210  double VertPosMean=0; //average position of all collection plane hits whose peaktime are in the window
211  double RMS=0; //rms of the collection plane hit positions
212  unsigned int CollectionHits=0; //total no of collection plane hits in the time window
213  for (unsigned int z=0; z < hitsZ.size(); ++z){
214  double PeakTimeMinusZ=hitsZ[z]->PeakTimePlusRMS(-1.);
215  //double PeakTimeZ=hitsZ[z]->PeakTime();
216  double ChargePerNumZ=hitsZ[z]->Integral()/hitsZ[z]->Multiplicity();
217  double rChargeZ=hitsUV[uv0]->Integral()/hitsZ[z]->Integral();
218  if (hitsUV[uv0]->Multiplicity()>0 && hitsZ[z]->Integral()>0) rChargeZ=ChargePerNum0/ChargePerNumZ;
219  double vert[3]={0.0,0.0,0.0};
220  std::vector<geo::WireID> wires = geo->ChannelToWire(hitsZ[z]->Channel());
221  geo->Cryostat(wires[0].Cryostat).TPC(wires[0].TPC).Plane(2).Wire(wires[0].Wire).GetCenter(vert);
222  if (fabs(fabs(PeakTimeMinusUV-PeakTimeMinusZ)-MinPeakTime)<10 &&
223  //fabs(PeakTimeUV-PeakTimeZ)<20 &&
224  (rChargeZ<rmax && rChargeZ>(1/rmax)) &&
225  wireiduv0[wireid0].TPC==wires[0].TPC ){
226  CollectionHits++;
227  VertPosMean+=vert[2];
228  RMS+=(VertPos-vert[2])*(VertPos-vert[2]);
229  }
230  }
231  VertPosMean=VertPosMean/double(CollectionHits); //average position
232  RMS=RMS/double(CollectionHits); //rms of the collection plane hits
233 
234 
235  //loop over the wireid of the induction plane hit whose time is closest to the hit in the first loop
236  //and find intersection and find a wire whose returned position is closest to the position of the collection
237  //plane hit found in the loop above.
238  for (unsigned int wireid1=0; wireid1 < wireiduv1.size(); ++wireid1){
239  if (abs(int(wireid0)-int(wireid1))<3 && wireiduv0[wireid0].TPC==wireiduv1[wireid1].TPC){
240  geo->IntersectionPoint(wireiduv0[wireid0].Wire, wireiduv1[wireid1].Wire,
241  wireiduv0[wireid0].Plane, wireiduv1[wireid1].Plane,
242  wireiduv0[wireid0].Cryostat, wireiduv0[wireid0].TPC,
243  CandPosition[1], CandPosition[2]);
244  //double VertBetUVandZ=fabs(VertPosMean-CandPosition[2]);
245  double VertBetUVandZ=fabs(VertPos-CandPosition[2]);
246  if (CandPosition[1]>=ybound[0] && CandPosition[1]<=ybound[1] && CandPosition[2]>=zbound[0] && CandPosition[2]<=zbound[1]){
247  if (VertBetUVandZ<MinVertBetUVandZ){
248  MinVertBetUVandZ=VertBetUVandZ;
249  FinalWireID=wireid0;
250  //FinalWire=wireiduv0[wireid0].Wire;
251  FinalPosition[1]=CandPosition[1];
252  FinalPosition[2]=CandPosition[2];
253  }
254  }
255  }
256  }
257  }
258 
259  //read simulated position of the induction plane hit using backtracker
260  // trj: skip this as it is part of the analysis and not the reco alg
261  //std::vector<double> SimPosition;
262  //unsigned int SimTPC, SimCstat, SimWire, SimWireID;
263  //SimPosition=bt->HitToXYZ(hitsUV[uv0]);
264  //if (SimPosition[2]!=99999){
265  // double xyzpos[3];
266  //xyzpos[0]=SimPosition[0];
267  //xyzpos[1]=SimPosition[1];
268  //xyzpos[2]=SimPosition[2];
269  // geo->PositionToTPC(xyzpos, SimTPC, SimCstat);
270  //SimWire=int(0.5+geo->WireCoordinate(xyzpos[1], xyzpos[2], hitsUV[uv0]->WireID().Plane, SimTPC, SimCstat));
271  // for (unsigned int wireid0=0; wireid0 < wireiduv0.size(); ++wireid0){
272  // if (abs(int(wireiduv0[wireid0].Wire)-int(SimWire))<2) SimWireID=wireid0;
273  // }
274  //}
275 
276  //save timeoffset corrected peaktime
277  double PeakTime=hitsUV[uv0]->PeakTime();
278  if (hitsUV[uv0]->View()==geo::kU) PeakTime+=timeoffset[0];
279  if (hitsUV[uv0]->View()==geo::kV) PeakTime+=timeoffset[1];
280 
281  //save hit information in a struct of a vector including positions returned in the disambiguation module
282  //and simulated hit information
283  h.Channel=hitsUV[uv0]->Channel();
284  h.StartTick=0;
285  h.EndTick=0;
286  h.PeakTime=hitsUV[uv0]->PeakTime();
287  h.SigmaPeakTime=hitsUV[uv0]->SigmaPeakTime();
288  h.RMS=hitsUV[uv0]->RMS();
289  h.PeakAmplitude=hitsUV[uv0]->PeakAmplitude();
290  h.SigmaPeakAmplitude=hitsUV[uv0]->SigmaPeakAmplitude();
291  h.SummedADC=hitsUV[uv0]->SummedADC();
292  h.Integral=hitsUV[uv0]->Integral();
293  h.SigmaIntegral=hitsUV[uv0]->SigmaIntegral();
294  h.Multiplicity=hitsUV[uv0]->Multiplicity();
295  h.LocalIndex=hitsUV[uv0]->LocalIndex();
296  h.GoodnessOfFit=hitsUV[uv0]->GoodnessOfFit();
297  h.DegreesOfFreedom=hitsUV[uv0]->DegreesOfFreedom();
298  h.View=hitsUV[uv0]->View();
299  h.SignalType=hitsUV[uv0]->SignalType();
300  h.WireID=wireiduv0[FinalWireID];
301  h.FinalYPos=FinalPosition[1];
302  h.FinalZPos=FinalPosition[2];
303  h.TPC=hitsUV[uv0]->WireID().TPC;
304  h.DisambigWireID=FinalWireID;
305  //h.SimYPos=SimPosition[1];
306  //h.SimZPos=SimPosition[2];
307  //h.SimWireID=SimWireID;
308  InductionAll.push_back(h);
309 
310  }//end of induction plane hit loop
311 
312 
313  //the disambiguation algorithm above returns a vector whose returned intersection position of two
314  //induction plane hits and returns a wireid of a wire whose position is closest to the returned
315  //position. the function below is to correct the position of hits whose side is correct disambiguated
316  //but the position is off. in the function, the detector volumne is going to be divided in xz position
317  //space and yposition of induction plane hits in the position space is going to be saved in a histogram,
318  //find a bin whose bincontents is highest, calculate the mean position as fitting the histogram around
319  //the bin, find a bin whose distance to the mean position found is about 100 and correct the position
320  //and wireid of hits in the bin.
321  int BinPerPos=4;
322  unsigned int CellsPerTPC=BinPerPos*BinPerPos;//cells in xz position space
323  double TotalXbin=(xbound[1]-xbound[0])/double(BinPerPos);
324  double TotalZbin=(zbound[1]-zbound[0])/double(BinPerPos);
325  //define histograms to save y position returned in the disambiguation algorithm above for all cells in
326  //xz position space
327  std::vector<TH1F *> YPosXZ(CellsPerTPC,0);
328  for (unsigned int cell=0; cell<YPosXZ.size(); cell++){
329  YPosXZ[cell]=new TH1F(Form("cell%d",cell),Form("cell%d",cell), 40, ybound[0], ybound[1]);
330  }
331 
332  //loop over all induction plane hits and save the returned yposition in a histogram
333  for (unsigned int uv0=0; uv0 < InductionAll.size(); ++uv0){
334  int XbinForHit=int( (InductionAll[uv0].PeakTime-xbound[0])/TotalXbin );
335  int ZbinForHit=int( (InductionAll[uv0].FinalZPos-zbound[0])/TotalZbin );
336  if (InductionAll[uv0].PeakTime>xbound[1]) XbinForHit=BinPerPos-1;
337  if (InductionAll[uv0].FinalZPos<zbound[0]){
338  ZbinForHit=0;
339  } else if (InductionAll[uv0].FinalZPos>zbound[1]){
340  ZbinForHit=BinPerPos-1;
341  }
342  int BinForHit=XbinForHit*BinPerPos+ZbinForHit;
343 
344  YPosXZ[BinForHit]->Fill(InductionAll[uv0].FinalYPos);
345  }//end of induction plane hit loop
346 
347  //define variables to save minimum and maximum range of the bin whose content is highest,
348  //a variable to indicate whether the position-correction function should be applied or not
349  //and variables to save minimum and maximum range of the bin whose mean position to the
350  //position of the highestcontent bin is 100cm away.
351  std::vector<double> FirstPeakXmin(CellsPerTPC,0);
352  std::vector<double> FirstPeakXmax(CellsPerTPC,0);
353  std::vector<double> ShouldCorrect(CellsPerTPC,0);
354  std::vector<double> SecondPeakXmin(CellsPerTPC,0);
355  std::vector<double> SecondPeakXmax(CellsPerTPC,0);
356 
357  //loop over all histograms
358  for (unsigned int cell=0; cell<CellsPerTPC; cell++){
359  //define variables to save highest bincontent and the bin index
360  double HitsInFirstPeak=0;
361  int FirstPeakBin=0;
362  for (unsigned int bin=0; bin<40; ++bin){
363  if (HitsInFirstPeak<YPosXZ[cell]->GetBinContent(bin+1)){
364  HitsInFirstPeak=YPosXZ[cell]->GetBinContent(bin+1);
365  FirstPeakBin=bin;
366  }
367  }
368 
369  //define the minimum and maximum range around the highest content bin
370  double PeakBinXmin=double(FirstPeakBin)*((ybound[1]-ybound[0])/40.0)+ybound[0]-10;
371  double PeakBinXmax=double(FirstPeakBin)*((ybound[1]-ybound[0])/40.0)+ybound[0]+10;
372 
373  //define a guassian histogram to fit the yposition histogram in the range calculated above
374  TF1 *YPosFitFirstPeak = new TF1("YPosFitFirstPeak","gaus", PeakBinXmin, PeakBinXmax);
375  YPosFitFirstPeak->SetParameter(0, 100);
376  YPosFitFirstPeak->SetParameter(1, PeakBinXmin+10);
377  YPosFitFirstPeak->SetParameter(2, 0.1);
378  //YPosFitFirstPeak->SetLineColor(4);
379  YPosXZ[cell]->Fit("YPosFitFirstPeak", "Q", "", PeakBinXmin-10, PeakBinXmax+10);
380  //save interal, mean and rms in variables
381  double FirstPeakBinAll=YPosFitFirstPeak->GetParameter(0);
382  double FirstPeakBinMean=YPosFitFirstPeak->GetParameter(1);
383  double FirstPeakBinRMS=YPosFitFirstPeak->GetParameter(2);
384 
385  FirstPeakXmin[cell]=FirstPeakBinMean-fabs(FirstPeakBinRMS)*2;
386  FirstPeakXmax[cell]=FirstPeakBinMean+fabs(FirstPeakBinRMS)*2;
387 
388  //define the position range where most incorrectly disambiguated hits are residing.
389  //depending on whether the mean position returned in above guassian fit is above or below
390  //the midpoint of the detector volume in ydirection, the position range could be 100cm below
391  //or above the mean position returned.
392  if (FirstPeakBinMean>(ybound[0]+ybound[1])*0.5){
393  SecondPeakXmin[cell]=FirstPeakXmin[cell]-100;
394  SecondPeakXmax[cell]=FirstPeakXmax[cell]-100;
395  } else if (FirstPeakBinMean<(ybound[0]+ybound[1])*0.5){
396  SecondPeakXmin[cell]=FirstPeakXmin[cell]+100;
397  SecondPeakXmax[cell]=FirstPeakXmax[cell]+100;
398  }
399 
400  //define a guassian histogram to fit the yposition histogram in the range where most incorrectly
401  //disambiguated hits are residing
402  TF1 *YPosFitSecondPeak = new TF1("YPosFitSecondPeak","gaus", SecondPeakXmin[cell], SecondPeakXmax[cell]);
403  YPosFitSecondPeak->SetParameter(0, 100);
404  YPosFitSecondPeak->SetParameter(1, SecondPeakXmin[cell]+10);
405  YPosFitSecondPeak->SetParameter(2, 0.1);
406  //YPosFitSecondPeak->SetLineColor(4);
407  YPosXZ[cell]->Fit("YPosFitSecondPeak", "Q", "", SecondPeakXmin[cell]-20, SecondPeakXmax[cell]+20);
408  //define variables to save integral, mean and rms of the guassian fit result.
409  double SecondPeakBinAll=YPosFitSecondPeak->GetParameter(0);
410  double SecondPeakBinMean=YPosFitSecondPeak->GetParameter(1);
411  double SecondPeakBinRMS=YPosFitSecondPeak->GetParameter(2);
412 
413  SecondPeakXmin[cell]=SecondPeakBinMean-fabs(SecondPeakBinRMS)*4;
414  SecondPeakXmax[cell]=SecondPeakBinMean+fabs(SecondPeakBinRMS)*4;
415 
416  //we assumed that the position-correction function works when more than half of the induction plane
417  //hits in the xz position cell is correctly disambiguated, meaning that most hits in the bin range whose
418  //bincontents is highest is correctly disambiguated. in this case, the hits in the highest bincontent
419  //ranges are in the firstpeak range and the incorrectly disambiguated hits in the secondpeak range. we
420  //only correct the position of the hits in the secondpeak range if the hits in the firstpeak is more than 5.
421  ShouldCorrect[cell]=0;
422  if (FirstPeakBinAll>SecondPeakBinAll &&
423  HitsInFirstPeak>5 &&
424  SecondPeakBinAll>0 && SecondPeakBinAll!=100 &&
425  SecondPeakXmin[cell]>ybound[0] && SecondPeakXmax[cell]<ybound[1]){
426  ShouldCorrect[cell]=1;
427  }
428 
429  //delete histograms used
430  YPosFitFirstPeak->Delete();
431  YPosFitSecondPeak->Delete();
432  YPosXZ[cell]->Delete();
433 
434  }//end of histogram loop
435 
436 
437  //now we know the hits whose position and wireid to be corrected so loop over the induction plane hits
438  //again, find the index of the hit in the xz position space and correct the position and wireid if the position
439  //returned in the disambiguation algorithm above is in the range where most incorrectly disambiguated hits could be residing.
440  for (unsigned int uv0=0; uv0 < InductionAll.size(); ++uv0){
441  int XbinForHit=int( (InductionAll[uv0].PeakTime-xbound[0])/TotalXbin );
442  int ZbinForHit=int( (InductionAll[uv0].FinalZPos-zbound[0])/TotalZbin );
443  if (InductionAll[uv0].PeakTime>xbound[1]) XbinForHit=BinPerPos-1;
444  if (InductionAll[uv0].FinalZPos<zbound[0]){
445  ZbinForHit=0;
446  } else if (InductionAll[uv0].FinalZPos>zbound[1]){
447  ZbinForHit=BinPerPos-1;
448  }
449  int BinForHit=XbinForHit*BinPerPos+ZbinForHit;//bin index for the hit
450 
451 
452  //shift the position and wireidreturned in the disambiguation algorithm
453  std::vector<geo::WireID> wireiduv0 = geo->ChannelToWire(InductionAll[uv0].Channel);
454  if (ShouldCorrect[BinForHit]!=0){
455  if (InductionAll[uv0].FinalYPos>SecondPeakXmin[BinForHit] &&
456  InductionAll[uv0].FinalYPos<SecondPeakXmax[BinForHit]){
457  InductionAll[uv0].FinalYPos+=InductionAll[uv0].FinalYPos+
458  (FirstPeakXmin[BinForHit]+FirstPeakXmax[BinForHit])*0.5-(SecondPeakXmin[BinForHit]+SecondPeakXmax[BinForHit])*0.5;
459 
460  if (InductionAll[uv0].DisambigWireID>1){
461  InductionAll[uv0].DisambigWireID=InductionAll[uv0].DisambigWireID-2;
462  } else if (InductionAll[uv0].DisambigWireID<2 && InductionAll[uv0].DisambigWireID+2<wireiduv0.size()){
463  InductionAll[uv0].DisambigWireID=InductionAll[uv0].DisambigWireID+2;
464  }
465  }
466  }
467 
468  //save the hit information in a hit object
469  recob::Hit hit(InductionAll[uv0].Channel,
470  InductionAll[uv0].StartTick,
471  InductionAll[uv0].EndTick,
472  InductionAll[uv0].PeakTime,
473  InductionAll[uv0].SigmaPeakTime,
474  InductionAll[uv0].RMS,
475  InductionAll[uv0].PeakAmplitude,
476  InductionAll[uv0].SigmaPeakAmplitude,
477  InductionAll[uv0].SummedADC,
478  InductionAll[uv0].Integral,
479  InductionAll[uv0].SigmaIntegral,
480  InductionAll[uv0].Multiplicity,
481  InductionAll[uv0].LocalIndex,
482  InductionAll[uv0].GoodnessOfFit,
483  InductionAll[uv0].DegreesOfFreedom,
484  InductionAll[uv0].View,
485  InductionAll[uv0].SignalType,
486  wireiduv0[InductionAll[uv0].DisambigWireID]);
487 
488  DisambiguatedHits.push_back(hit);
489 
490  //if (InductionAll[uv0].SimYPos!=99999){
491  //HitsTotal+=InductionAll[uv0].Integral;
492  //if (InductionAll[uv0].DisambigWireID==InductionAll[uv0].SimWireID) HitsCorrect+=InductionAll[uv0].Integral;
493  //}
494 
495  }//end of induction plane hit loop
496 
497  }
498  }
499 
500  //cout << "testing dune disambig" << "efficiency:" << double(HitsCorrect)*100.0/double(HitsTotal) << endl;
501 
502 
503 }
504 
505 } //end namespace apa
double PeakTime
unsigned int EndTick
double FinalZPos
unsigned int DisambigWireID
Encapsulate the construction of a single cyostat.
AdcChannelData::View View
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
Planes which measure V.
Definition: geo_types.h:130
double SummedADC
unsigned int Channel
struct vector vector
double SigmaPeakTime
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
STL namespace.
Planes which measure Z direction.
Definition: geo_types.h:132
int DegreesOfFreedom
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
double SigmaIntegral
double FinalYPos
unsigned int LocalIndex
double Length() const
Length is associated with z coordinate [cm].
Definition: TPCGeo.h:115
double Integral
T abs(T value)
Planes which measure U.
Definition: geo_types.h:129
unsigned int StartTick
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
virtual void reconfigure(fhicl::ParameterSet const &pset)
double GoodnessOfFit
geo::View_t View
geo::SigType_t SignalType
T get(std::string const &key) const
Definition: ParameterSet.h:271
p
Definition: test.py:223
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
double PeakAmplitude
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
bool IntersectionPoint(geo::WireID const &wid1, geo::WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
Encapsulate the geometry of a wire.
double HalfHeight() const
Height is associated with y coordinate [cm].
Definition: TPCGeo.h:111
Declaration of signal hit object.
Encapsulate the construction of a single detector plane.
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
unsigned int TPC
QTextStream & bin(QTextStream &s)
geo::WireID WireID
double SigmaPeakAmplitude
Declaration of basic channel signal object.
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
recob::tracking::Plane Plane
Definition: TrackState.h:17
LArSoft geometry interface.
Definition: ChannelGeo.h:16
void LocalToWorld(const double *tpc, double *world) const
Transform point from local TPC frame to world frame.
Definition: TPCGeo.h:563
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
unsigned int Multiplicity
Encapsulate the construction of a single detector plane.
double RMS