SignalToNoise_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SignalToNoise
3 // Module Type: analyzer
4 // File: SignalToNoise_module.cc
5 //
6 // Generated at Sat Sep 10 22:01:34 2016 by Tingjun Yang using artmod
7 // from cetpkgsupport v1_10_02.
8 ////////////////////////////////////////////////////////////////////////
9 
16 #include "art_root_io/TFileService.h"
18 #include "fhiclcpp/ParameterSet.h"
34 #include "TH1D.h"
35 #include "TH2D.h"
36 #include "TGraph.h"
37 #include "TTree.h"
38 
39 #include <vector>
40 
41 namespace dune {
42  class SignalToNoise;
43 }
44 
46 public:
47  explicit SignalToNoise(fhicl::ParameterSet const & p);
48  // The destructor generated by the compiler is fine for classes
49  // without bare pointers or other resource use.
50 
51  // Plugins should not be copied or assigned.
52  SignalToNoise(SignalToNoise const &) = delete;
53  SignalToNoise(SignalToNoise &&) = delete;
54  SignalToNoise & operator = (SignalToNoise const &) = delete;
56 
57  // Required functions.
58  void analyze(art::Event const & e) override;
59 
60  // Selected optional functions.
61  void beginJob() override;
62 
63 private:
64 
65  // Declare member data here.
70 
71  bool fIsData;
72 
74 
75  double cx[93], cy[93], cz[93];
76 
77  TH1D *dcos;
78  TH2D *hyz, *hxz;
79  TH1D *hsig[8][3][4];
80  TH1D *hbkg[8][3][4];
81  TH1D *hstb[8][3][4];
82  TH1D *hdx;
83  TH2D *hphx;
84  TH1D *hdt;
85 
86  TTree *ftree;
87  int run;
88  int event;
89  int tpc;
90  int wire;
91  int trackid;
92  double dqdx;
93  double x;
94 };
95 
96 
98  :
99  EDAnalyzer(p),
100  fTrackModuleLabel(p.get< art::InputTag >("TrackModuleLabel")),
101  fExternalCounterModuleLabel(p.get< art::InputTag >("ExternalCounterModuleLabel")),
102  fRawDigitModuleLabel(p.get< art::InputTag >("RawDigitModuleLabel")),
103  fCalDataModuleLabel(p.get< art::InputTag >("CalDataModuleLabel")),
104  fIsData(p.get<bool>("IsData",true))
105 {
106 }
107 
109 {
110  // Implementation of required member function here.
111  run = e.run();
112  event = e.id().event();
113 
114  // * tracks
115  std::vector<art::Ptr<recob::Track> > tracklist;
116  auto trackListHandle = e.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
117  if (trackListHandle)
118  art::fill_ptr_vector(tracklist, trackListHandle);
119 
120  // * External Counters
121  std::vector<art::Ptr<raw::ExternalTrigger> > countlist;
122  auto countListHandle = e.getHandle< std::vector<raw::ExternalTrigger> >(fExternalCounterModuleLabel);
123  if (countListHandle)
124  art::fill_ptr_vector(countlist, countListHandle);
125 
126  // * Raw Digits
127  std::vector<art::Ptr<raw::RawDigit> > rawlist;
128  auto rawListHandle = e.getHandle<std::vector<raw::RawDigit> >(fRawDigitModuleLabel);
129  if (rawListHandle)
130  art::fill_ptr_vector(rawlist, rawListHandle);
131 
132  // * Noise removed
133  std::vector<art::Ptr<recob::Wire> > wirelist;
134  auto wireListHandle = e.getHandle<std::vector<recob::Wire> >(fCalDataModuleLabel);
135  if (wireListHandle)
136  art::fill_ptr_vector(wirelist, wireListHandle);
137 
138  // * associations
139  art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, e, fTrackModuleLabel);
140 
143 
144  std::vector<int> c1; //counter index
145  std::vector<int> c2;
146  std::vector<int> ci1; //trigger id
147  std::vector<int> ci2;
148 
149  for (size_t i = 0; i<countlist.size(); ++i){
150  for (size_t j = 0; j<countlist.size(); ++j){
151  if (i==j) continue;
152  if (std::abs(countlist[i]->GetTrigTime()-countlist[j]->GetTrigTime()) < 2){
153  // for c2: GetTrigID() is an unsigned int and always >= 0
154  //if ( ( countlist[i]->GetTrigID() >= 6 && countlist[i]->GetTrigID() <= 15 && countlist[j]->GetTrigID() >= 28 && countlist[j]->GetTrigID() <= 37 ) // East Lower, West Upper
155  // || ( countlist[i]->GetTrigID() >= 0 && countlist[i]->GetTrigID() <= 5 && countlist[j]->GetTrigID() >= 22 && countlist[j]->GetTrigID() <= 27 ) // South Lower, North Upper
156  // || ( countlist[i]->GetTrigID() >= 16 && countlist[i]->GetTrigID() <= 21 && countlist[j]->GetTrigID() >= 38 && countlist[j]->GetTrigID() <= 43 ) // North Lower, South Upper
157  if ( ( countlist[i]->GetTrigID() >= 6 && countlist[i]->GetTrigID() <= 15 && countlist[j]->GetTrigID() >= 28 && countlist[j]->GetTrigID() <= 37 ) // East Lower, West Upper
158  || ( countlist[i]->GetTrigID() <= 5 && countlist[j]->GetTrigID() >= 22 && countlist[j]->GetTrigID() <= 27 ) // South Lower, North Upper
159  || ( countlist[i]->GetTrigID() >= 16 && countlist[i]->GetTrigID() <= 21 && countlist[j]->GetTrigID() >= 38 && countlist[j]->GetTrigID() <= 43 ) // North Lower, South Upper
160  ) {
161  // std::cout << "I have a match..."
162  // << "i " << i << ", ID " << countlist[i]->GetTrigID() << ", time " << countlist[i]->GetTrigTime() << "..."
163  // << "j " << j << ", ID " << countlist[j]->GetTrigID() << ", Time " << countlist[j]->GetTrigTime()
164  // << std::endl;
165  c1.push_back(countlist[i]->GetTrigID());
166  c2.push_back(countlist[j]->GetTrigID());
167  ci1.push_back(i);
168  ci2.push_back(j);
169  }
170  }
171  }
172  }
173  std::cout<<"counter array size = "<<c1.size()<<" tracklist size = "<<tracklist.size()<<std::endl;
174  if (c1.size()&&tracklist.size()){
175  double ccosx = (cx[c1[0]]-cx[c2[0]])/sqrt(pow(cx[c1[0]]-cx[c2[0]],2)+pow(cy[c1[0]]-cy[c2[0]],2)+pow(cz[c1[0]]-cz[c2[0]],2));
176  double ccosy = (cy[c1[0]]-cy[c2[0]])/sqrt(pow(cx[c1[0]]-cx[c2[0]],2)+pow(cy[c1[0]]-cy[c2[0]],2)+pow(cz[c1[0]]-cz[c2[0]],2));
177  double ccosz = (cz[c1[0]]-cz[c2[0]])/sqrt(pow(cx[c1[0]]-cx[c2[0]],2)+pow(cy[c1[0]]-cy[c2[0]],2)+pow(cz[c1[0]]-cz[c2[0]],2));
178 
179  recob::Track::Vector_t larStart;
180  recob::Track::Vector_t larEnd;
181 
182  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e);
183  for (size_t i = 0; i<tracklist.size(); ++i){
184  recob::Track::Point_t trackStart, trackEnd;
185  std::tie(trackStart, trackEnd) = tracklist[i]->Extent();
186  larStart = tracklist[i]->VertexDirection();
187  larEnd = tracklist[i]->EndDirection();
188 
189  double dc = std::abs(ccosx*larStart.X()+ccosy*larStart.Y()+ccosz*larStart.Z());
190  dcos->Fill(dc);
191  /*
192  if (i==4){
193  for (size_t j = 0; j<tracklist[i]->NPoints(); ++j){
194  if (tracklist[i]->HasValidPoint(j)){
195  std::cout<<j<<" "<<tracklist[i]->LocationAtPoint(j).X()
196  <<" "<<tracklist[i]->LocationAtPoint(j).Y()
197  <<" "<<tracklist[i]->LocationAtPoint(j).Z()
198  <<" "<<tracklist[i]->DirectionAtPoint(j).X()
199  <<" "<<tracklist[i]->DirectionAtPoint(j).Y()
200  <<" "<<tracklist[i]->DirectionAtPoint(j).Z()
201  <<std::endl;
202  }
203  }
204  }
205  */
206  if (dc>0.98){//found a match
207  double ctime = (countlist[ci1[0]]->GetTrigTime() + countlist[ci2[0]]->GetTrigTime())/(2*32.); //ticks
208  double x0 = trackStart.X() + (trackStart.X()>0?-1:1)*ctime*0.5*0.1085;
209  //double y0 = trackStart.Y();
210  double z0 = trackStart.Z();
211  double x1 = trackEnd.X() + (trackStart.X()>0?-1:1)*ctime*0.5*0.1085;
212  //double y1 = trackEnd.Y();
213  double z1 = trackEnd.Z();
214 
215  double dx1 = 1e10;
216  double dx2 = 1e10;
217  if (c1[0]>=6&&c1[0]<=15){
218  dx1 = std::abs(x0+(cz[c1[0]]-z0)*(x0-x1)/(z0-z1)-cx[c1[0]]);
219  }
220  if (c2[0]>=28&&c2[0]<=37){
221  dx2 = std::abs(x0+(cz[c2[0]]-z0)*(x0-x1)/(z0-z1)-cx[c2[0]]);
222  }
223  //std::cout<<dx1<<" "<<dx2<<std::endl;
224  hdx->Fill(x0+(cz[c1[0]]-z0)*(x0-x1)/(z0-z1)-cx[c1[0]]);
225  hdx->Fill(x0+(cz[c2[0]]-z0)*(x0-x1)/(z0-z1)-cx[c2[0]]);
226  if (dx1>40||dx2>40) continue;
227  for (size_t j = 0; j<tracklist[i]->NumberTrajectoryPoints(); ++j){
228  auto loc = tracklist[i]->LocationAtPoint(j);
229  hyz->Fill(loc.Z(),loc.Y());
230  if (loc.X()>0){
231  hxz->Fill(loc.Z(),loc.X()-ctime*0.5*0.1085);
232  }
233  else{
234  hxz->Fill(loc.Z(),loc.X()+ctime*0.5*0.1085);
235  }
236  }
237  //map wireid with hit pointer, index on track and x position
238  std::map<size_t,art::Ptr<recob::Hit> > hitmap;
239  std::map<size_t,int> indexmap;
240  std::map<size_t,double> xposmap;
241  if (fmthm.isValid()){//fmthm.isValid()
242  auto vhit = fmthm.at(i);
243  auto vmeta = fmthm.data(i);
244 
245  std::map<geo::PlaneID, std::vector<std::pair<int, double>>> planehits;
246 
247  for (size_t h = 0; h < vhit.size(); ++h){
248  auto loc = tracklist[i]->LocationAtPoint(vmeta[h]->Index());
249  if (vhit[h]->WireID().TPC==5&&vhit[h]->WireID().Plane==2){
250  hitmap[vhit[h]->WireID().Wire] = vhit[h];
251  indexmap[vhit[h]->WireID().Wire] = vmeta[h]->Index();
252  xposmap[vhit[h]->WireID().Wire] = loc.X()-ctime*0.5*0.1085;
253  }
254  planehits[vhit[h]->WireID().planeID()].push_back(std::make_pair(vhit[h]->WireID().Wire, vhit[h]->PeakTime()));
255  //now study the difference between hit time and prediction from counters
256  double xyzStart[3], xyzEnd[3]; //wire ends
257  geom->WireEndPoints(vhit[h]->WireID(), xyzStart, xyzEnd);
258  double xhit, yhit, zhit; //intersection between wire and counter line
259  //std::cout<<"about to calculate intersection"<<std::endl;
260  //std::cout<<xyzStart[1]<<" "<<xyzStart[2]<<" "<<xyzEnd[1]<<" "<<xyzEnd[2]<<" "<<cy[c1[0]]<<" "<<cz[c1[0]]<<" "<<cy[c2[0]]<<" "<<cz[c2[0]]<<std::endl;
261  if (geom->IntersectLines(xyzStart[1], xyzStart[2],
262  xyzEnd[1], xyzEnd[2],
263  cy[c1[0]], cz[c1[0]],
264  cy[c2[0]], cz[c2[0]],
265  yhit, zhit)){
266  xhit = cx[c1[0]] + (cx[c2[0]]-cx[c1[0]])/(cy[c2[0]]-cy[c1[0]])*(yhit-cy[c1[0]]);
267  hdt->Fill(vhit[h]->PeakTime() - detProp.ConvertXToTicks(xhit+(vhit[h]->WireID().TPC%2==1?1:-1)*ctime*0.5*0.1085, vhit[h]->WireID()));
268  }
269  }//loop over all associated hits
270 
271  //make sure there are enough hits on each plane before proceeding
272  for (auto &plhits : planehits){//loop over planes
273  int w0 = INT_MAX;
274  int w1 = INT_MIN;
275  std::vector<float> x, y, ey;
276  for (auto &hit : plhits.second){
277  if (hit.first < w0) w0 = hit.first;
278  if (hit.first > w1) w1 = hit.first;
279  x.push_back(hit.first);
280  y.push_back(hit.second);
281  ey.push_back(10); //arbitrary error on time
282  }
283  if (w0!=INT_MAX&&w1!=INT_MIN&&(w1-w0>5)&&plhits.second.size()>2){//there are at least 3 hits and the track is at least 5 wires long
284  float intcpt, slope, intcpterr, slopeerr, chidof;
285  fLinFitAlg.LinFit(x, y, ey, intcpt, slope, intcpterr, slopeerr, chidof);
286  //Loop over all wires in the range identified
287  std::map<raw::ChannelID_t, int> chused;
288  for (int w = w0; w<=w1; ++w){
289  geo::WireID wid(plhits.first, w);
290  auto channel = geom->PlaneWireToChannel(wid);
291  if (fCSP->IsBad(channel)) continue;
292  if (chused.find(channel)==chused.end()){ //channel not looked at before
293  double tick = intcpt + slope*w;
294  double x = detProp.ConvertTicksToX(tick, wid);//raw x
295  x += (wid.TPC%2==1?-1:1)*ctime*0.5*0.1085; //correct for t0
296  double xyzStart[3], xyzEnd[3]; //wire ends
297  geom->WireEndPoints(wid, xyzStart, xyzEnd);
298  x = std::abs(x-xyzStart[0]); //distance to wire plane
299  if (x<50){//hit is within 50 cm of wire plane
300  int besttime = -1;
301  for (size_t j = 0; j<rawlist.size(); ++j){
302  if (rawlist[j]->Channel() == channel){
303  double pedestal = rawlist[j]->GetPedestal();
304  int maxph = -1;
305  for (size_t k = 0; k<rawlist[j]->NADC(); ++k){
306  if (float(k)>=tick&&float(k)<=tick+20){
307  if (int(rawlist[j]->ADC(k)-pedestal)>maxph){
308  maxph = int(rawlist[j]->ADC(k)-pedestal);
309  besttime = k;
310  }
311  }
312  }
313  double mean=0;
314  double mean2=0;
315  int npts = 0;
316  for (size_t k = 0; k<rawlist[j]->NADC(); ++k){
317  if ((int(k)>besttime-200&&int(k)<besttime-100)||
318  (int(k)>besttime+100&&int(k)<besttime+200)){
319  mean += rawlist[j]->ADC(k)-pedestal;
320  mean2 += pow(rawlist[j]->ADC(k)-pedestal,2);
321  ++npts;
322  }
323  }
324  mean/=npts;
325  mean2/=npts;
326  double angleToVert = geom->WireAngleToVertical(geom->View(wid), wid.TPC, wid.Cryostat) - 0.5*::util::pi<>();
327  //std::cout<<vhit[h]->View()<<" "<<vhit[h]->WireID().TPC<<" "<<vhit[h]->WireID().Cryostat<<" "<<angleToVert<<std::endl;
328  const auto& dir = tracklist[i]->DirectionAtPoint(0);
329  double cosgamma = std::abs(std::sin(angleToVert)*dir.Y() + std::cos(angleToVert)*dir.Z());
330  //std::cout<<maxph*cosgamma<<" "<<sqrt(mean2-mean*mean)<<std::endl;
331  //std::cout<<vhit[h]->PeakTime()<<" "<<vhit[h]->PeakAmplitude()<<" "<<besttime<<" "<<maxph<<std::endl;
332  hsig[wid.TPC][wid.Plane][0]->Fill(maxph*cosgamma);
333  hbkg[wid.TPC][wid.Plane][0]->Fill(sqrt(mean2-mean*mean));
334  }
335  }//loop over rawlist
336  if (besttime<0||besttime>=15000) continue;
337  for (size_t j = 0; j<wirelist.size(); ++j){
338  if (wirelist[j]->Channel() == channel){
339  double maxph1 = wirelist[j]->Signal()[besttime];
340  double mean=0;
341  double mean2=0;
342  int npts = 0;
343  for (size_t k = 0; k<wirelist[j]->NSignal(); ++k){
344  if ((int(k)>besttime-200&&int(k)<besttime-100)||
345  (int(k)>besttime+100&&int(k)<besttime+200)){
346  mean += wirelist[j]->Signal()[k];
347  mean2 += pow(wirelist[j]->Signal()[k],2);
348  ++npts;
349  }
350  }
351  mean/=npts;
352  mean2/=npts;
353  double angleToVert = geom->WireAngleToVertical(geom->View(wid), wid.TPC, wid.Cryostat) - 0.5*::util::pi<>();
354  //std::cout<<vhit[h]->View()<<" "<<vhit[h]->WireID().TPC<<" "<<vhit[h]->WireID().Cryostat<<" "<<angleToVert<<std::endl;
355  const auto& dir = tracklist[i]->DirectionAtPoint(0);
356  double cosgamma = std::abs(std::sin(angleToVert)*dir.Y() + std::cos(angleToVert)*dir.Z());
357  //std::cout<<maxph*cosgamma<<" "<<sqrt(mean2-mean*mean)<<std::endl;
358  //std::cout<<vhit[h]->PeakTime()<<" "<<vhit[h]->PeakAmplitude()<<" "<<besttime<<" "<<maxph<<std::endl;
359  hsig[wid.TPC][wid.Plane][1]->Fill(maxph1*cosgamma);
360  hbkg[wid.TPC][wid.Plane][1]->Fill(sqrt(mean2-mean*mean));
361  }
362  }//loop over wirelist
363  }//within 50 cm of anode
364  chused[channel] = 1;
365  }//channel not used
366  }//loop over wires in range
367  }//find enough wires
368  }//loop over planes
369  }//fmthm is valid
370  if (hitmap.size()>=10){//found at least 10 hits in TPC5
371  for (size_t h = 0; h<geom->Nwires(2,5,0); ++h){//plane 2, tpc 5, cstat 0
372  if (fCSP->IsBad(geom->PlaneWireToChannel(2,h,5))) continue;
373  double pt = -1;
374  double xpos = -1;
375  if (hitmap.find(h)!=hitmap.end()){//found hit on track
376  pt = hitmap[h]->PeakTime();
377  xpos = xposmap[h];
378  }
379  else{//found hit time through extrapolation
380  std::vector<double> vw;
381  std::vector<double> vt;
382  std::vector<double> vx;
383  for (auto& hv : hitmap){
384  // for c2: fix ambiguous call to abs
385  //if (std::abs(hv.first-h)<=5){
386  if (util::absDiff(hv.first,h)<=5){
387  vw.push_back((hv.second)->WireID().Wire);
388  vt.push_back((hv.second)->PeakTime());
389  vx.push_back(xposmap[hv.first]);
390  }
391  }
392  if (vw.size()>=3){
393  TGraph *gr = new TGraph(vw.size(), &vw[0], &vt[0]);
394  pt = gr->Eval(h);
395  delete gr;
396  gr = new TGraph(vw.size(), &vw[0], &vx[0]);
397  xpos = gr->Eval(h);
398  delete gr;
399  }
400  }
401  if (pt>=0){
402  double maxph = -1e10;
403  for (size_t j = 0; j<rawlist.size(); ++j){
404  if (rawlist[j]->Channel() == geom->PlaneWireToChannel(2,h,5)){
405  double pedestal = rawlist[j]->GetPedestal();
406  for (size_t k = 0; k<rawlist[j]->NADC(); ++k){
407  if (float(k)>=pt&&float(k)<=pt+20){
408  if (rawlist[j]->ADC(k)-pedestal>maxph){
409  maxph = rawlist[j]->ADC(k)-pedestal;
410  }
411  }
412  }
413  }
414  }
415  double angleToVert = geom->WireAngleToVertical(hitmap.begin()->second->View(), hitmap.begin()->second->WireID().TPC, hitmap.begin()->second->WireID().Cryostat) - 0.5*::util::pi<>();
416  //std::cout<<vhit[h]->View()<<" "<<vhit[h]->WireID().TPC<<" "<<vhit[h]->WireID().Cryostat<<" "<<angleToVert<<std::endl;
417  //const TVector3& dir = tracklist[i]->DirectionAtPoint(indexmap.begin()->second);
418  const auto& dir = tracklist[i]->VertexDirection();
419  double cosgamma = std::abs(std::sin(angleToVert)*dir.Y() + std::cos(angleToVert)*dir.Z());
420  hphx->Fill(xpos,maxph*cosgamma/geom->WirePitch(hitmap.begin()->second->View()));
421  tpc = 5;
422  wire = h;
423  trackid = i;
424  x = xpos;
425  dqdx = maxph*cosgamma/geom->WirePitch(hitmap.begin()->second->View());
426  ftree->Fill();
427  //std::cout<<h<<" "<<pt<<" "<<xpos<<" "<<maxph*cosgamma/geom->WirePitch(hitmap.begin()->second->View())<<std::endl;
428  }
429  }
430  }
431  }//matched track
432  }
433  }
434 }
435 
437 {
439 
440  dcos = tfs->make<TH1D>("dcos",";cos#theta;",100,0,1.1);
441  dcos->Sumw2();
442 
443  hxz = tfs->make<TH2D>("hxz",";z (cm);x (cm)",1000,0,160,1000,-50,250);
444  hyz = tfs->make<TH2D>("hyz",";z (cm);y (cm)",1000,0,160,1000,-90,120);
445 
446  for (int i = 0; i<8; ++i){
447  for (int j = 0; j<3; ++j){
448  for (int k = 0; k<4; ++k){
449  hsig[i][j][k]= tfs->make<TH1D>(Form("hsig_%d_%d_%d",i,j,k),Form("TPC=%d, Plane=%d, %d",i,j,k), 200,0,200);
450  hbkg[i][j][k]= tfs->make<TH1D>(Form("hbkg_%d_%d_%d",i,j,k),Form("TPC=%d, Plane=%d, %d",i,j,k), 100,0,100);
451  hstb[i][j][k]= tfs->make<TH1D>(Form("hstb_%d_%d_%d",i,j,k),Form("TPC=%d, Plane=%d, %d",i,j,k), 100,0,100);
452  }
453  }
454  }
455 
456  hdx = tfs->make<TH1D>("hdx",";#Delta x (cm)",100,-100,100);
457  hdx->Sumw2();
458  hphx = tfs->make<TH2D>("hphx",";x (cm); dQ/dx (ADC/cm)",50,0,250,100,0,600);
459 
460  hdt = tfs->make<TH1D>("hdt",";#Delta t (ticks);",1000,-1000,1000);
461 
462  ftree = tfs->make<TTree>("tree","a tree for sn analysis");
463  ftree->Branch("run", &run, "run/I");
464  ftree->Branch("event", &event, "event/I");
465  ftree->Branch("tpc", &tpc, "tpc/I");
466  ftree->Branch("wire", &wire, "wire/I");
467  ftree->Branch("trackid", &trackid, "trackid/I");
468  ftree->Branch("dqdx", &dqdx, "dqdx/D");
469  ftree->Branch("x", &x, "x/D");
470 
471  if (fIsData){
472  std::ifstream in;
473  //in.open("/dune/app/users/mthiesse/olddev/CounterZOffset/work/counterInformation.txt");
474  in.open("counterInformation.txt");
475  char line[1024];
476  while(1){
477  in.getline(line,1024);
478  if (!in.good()) break;
479  int i;
480  float x,y,z;
481  sscanf(line,"%d %f %f %f",&i,&x,&y,&z);
482  //std::cout<<i<<" "<<x<<" "<<y<<" "<<z<<std::endl;
483  cx[i] = x;
484  cy[i] = y;
485  cz[i] = z;
486  }
487  in.close();
488  }
489  else{
491  for (size_t i = 0; i<geom->NAuxDets(); ++i){
492  auto& auxdet = geom->AuxDet(i);
493  //std::cout<<i<<" "<<auxdet.GetCenter().X()<<" "<<auxdet.GetCenter().Y()<<" "<<auxdet.GetCenter().Z()<<std::endl;
494  cx[i] = auxdet.GetCenter().X();
495  cy[i] = auxdet.GetCenter().Y();
496  cz[i] = auxdet.GetCenter().Z();
497  }
498  }
499 }
500 
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
Functions to help with numbers.
art::InputTag fTrackModuleLabel
TH3F * xpos
Definition: doAna.cpp:29
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
art::InputTag fExternalCounterModuleLabel
constexpr T pow(T x)
Definition: pow.h:72
Class to keep data related to recob::Hit associated with recob::Track.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
uint8_t channel
Definition: CRTFragment.hh:201
string dir
unsigned int Index
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: types.h:32
void analyze(art::Event const &e) override
trkf::LinFitAlg fLinFitAlg
unsigned int Nwires(unsigned int p, unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wires in the specified plane.
art framework interface to geometry description
tracking::Vector_t Vector_t
Definition: Track.h:54
void LinFit(std::vector< float > &x, std::vector< float > &y, std::vector< float > &ey2, float &Intercept, float &Slope, float &InterceptError, float &SlopeError, float &ChiDOF) const
Definition: LinFitAlg.cxx:17
T abs(T value)
const double e
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
p
Definition: test.py:223
art::InputTag fCalDataModuleLabel
Class providing information about the quality of channels.
RunNumber_t run() const
Definition: DataViewImpl.cc:71
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
SignalToNoise(fhicl::ParameterSet const &p)
bool IntersectLines(double A_start_x, double A_start_y, double A_end_x, double A_end_y, double B_start_x, double B_start_y, double B_end_x, double B_end_y, double &x, double &y) const
Computes the intersection between two lines on a plane.
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
AuxDetGeo const & AuxDet(unsigned int const ad=0) const
Returns the specified auxiliary detector.
Detector simulation of raw signals on wires.
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
Definition: NumericUtils.h:43
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Declaration of signal hit object.
void line(double t, double *p, double &x, double &y, double &z)
SignalToNoise & operator=(SignalToNoise const &)=delete
tracking::Point_t Point_t
Definition: Track.h:53
Provides recob::Track data product.
Interface for experiment-specific channel quality info provider.
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
art::InputTag fRawDigitModuleLabel
EventNumber_t event() const
Definition: EventID.h:116
Declaration of basic channel signal object.
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
Interface for experiment-specific service for channel quality info.
recob::tracking::Plane Plane
Definition: TrackState.h:17
int bool
Definition: qglobal.h:345
Utility functions to extract information from recob::Track
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
double WireAngleToVertical(geo::View_t view, geo::TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
EventID id() const
Definition: Event.cc:34
QTextStream & endl(QTextStream &s)
unsigned int NAuxDets() const
Returns the number of auxiliary detectors.
Event finding and building.