Public Member Functions | Private Attributes | List of all members
dune::SignalToNoise Class Reference
Inheritance diagram for dune::SignalToNoise:
art::EDAnalyzer art::detail::Analyzer art::detail::LegacyModule art::Observer art::ModuleBase

Public Member Functions

 SignalToNoise (fhicl::ParameterSet const &p)
 
 SignalToNoise (SignalToNoise const &)=delete
 
 SignalToNoise (SignalToNoise &&)=delete
 
SignalToNoise & operator= (SignalToNoise const &)=delete
 
SignalToNoise & operator= (SignalToNoise &&)=delete
 
void analyze (art::Event const &e) override
 
void beginJob () override
 
- Public Member Functions inherited from art::EDAnalyzer
 EDAnalyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDAnalyzer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Analyzer
virtual ~Analyzer () noexcept
 
 Analyzer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 Analyzer (Table< Config > const &config)
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Observer
 ~Observer () noexcept
 
 Observer (Observer const &)=delete
 
 Observer (Observer &&)=delete
 
Observer & operator= (Observer const &)=delete
 
Observer & operator= (Observer &&)=delete
 
void registerProducts (ProductDescriptions &, ModuleDescription const &)
 
void fillDescriptions (ModuleDescription const &)
 
fhicl::ParameterSetID selectorConfig () const
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Attributes

art::InputTag fTrackModuleLabel
 
art::InputTag fExternalCounterModuleLabel
 
art::InputTag fRawDigitModuleLabel
 
art::InputTag fCalDataModuleLabel
 
bool fIsData
 
trkf::LinFitAlg fLinFitAlg
 
double cx [93]
 
double cy [93]
 
double cz [93]
 
TH1D * dcos
 
TH2D * hyz
 
TH2D * hxz
 
TH1D * hsig [8][3][4]
 
TH1D * hbkg [8][3][4]
 
TH1D * hstb [8][3][4]
 
TH1D * hdx
 
TH2D * hphx
 
TH1D * hdt
 
TTree * ftree
 
int run
 
int event
 
int tpc
 
int wire
 
int trackid
 
double dqdx
 
double x
 

Additional Inherited Members

- Public Types inherited from art::EDAnalyzer
using WorkerType = WorkerT< EDAnalyzer >
 
using ModuleType = EDAnalyzer
 
- Protected Member Functions inherited from art::Observer
std::string const & processName () const
 
bool wantAllEvents () const noexcept
 
bool wantEvent (ScheduleID id, Event const &e) const
 
Handle< TriggerResults > getTriggerResults (Event const &e) const
 
 Observer (fhicl::ParameterSet const &config)
 
 Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 45 of file SignalToNoise_module.cc.

Constructor & Destructor Documentation

dune::SignalToNoise::SignalToNoise ( fhicl::ParameterSet const &  p)
explicit

Definition at line 97 of file SignalToNoise_module.cc.

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 }
dune::SignalToNoise::fTrackModuleLabel
art::InputTag fTrackModuleLabel
Definition: SignalToNoise_module.cc:66
dune::SignalToNoise::fExternalCounterModuleLabel
art::InputTag fExternalCounterModuleLabel
Definition: SignalToNoise_module.cc:67
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
test.p
p
Definition: test.py:223
dune::SignalToNoise::fCalDataModuleLabel
art::InputTag fCalDataModuleLabel
Definition: SignalToNoise_module.cc:69
dune::SignalToNoise::fRawDigitModuleLabel
art::InputTag fRawDigitModuleLabel
Definition: SignalToNoise_module.cc:68
dune::SignalToNoise::SignalToNoise ( SignalToNoise const &  )
delete
dune::SignalToNoise::SignalToNoise ( SignalToNoise &&  )
delete

Member Function Documentation

void dune::SignalToNoise::analyze ( art::Event const &  e)
overridevirtual

Implements art::EDAnalyzer.

Definition at line 108 of file SignalToNoise_module.cc.

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 
141  art::ServiceHandle<geo::Geometry> geom;
142  const lariov::ChannelStatusProvider* fCSP = &art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
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 }
lariov::ChannelStatusProvider::IsBad
virtual bool IsBad(raw::ChannelID_t channel) const =0
Returns whether the specified channel is bad in the current run.
dune::SignalToNoise::fTrackModuleLabel
art::InputTag fTrackModuleLabel
Definition: SignalToNoise_module.cc:66
xpos
TH3F * xpos
Definition: doAna.cpp:29
reco_momentum_tuples.float
float
Definition: reco_momentum_tuples.py:12
dune::SignalToNoise::fExternalCounterModuleLabel
art::InputTag fExternalCounterModuleLabel
Definition: SignalToNoise_module.cc:67
cet::pow
constexpr T pow(T x)
Definition: pow.h:72
channel
uint8_t channel
Definition: CRTFragment.hh:201
dir
string dir
Definition: INukeNucleonCorr.cxx:61
Index
unsigned int Index
Definition: sspmapmaker_v1.c:14
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:824
Signal
Definition: types.h:32
dune::SignalToNoise::fLinFitAlg
trkf::LinFitAlg fLinFitAlg
Definition: SignalToNoise_module.cc:73
geo::GeometryCore::Nwires
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.
Definition: GeometryCore.h:3342
recob::Track::Vector_t
tracking::Vector_t Vector_t
Definition: Track.h:54
trkf::LinFitAlg::LinFit
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
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
e
const double e
Definition: gUpMuFluxGen.cxx:165
geo::fhicl::WireID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Definition: geo_types_fhicl.h:407
util::quantities::tick
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
dune::SignalToNoise::fCalDataModuleLabel
art::InputTag fCalDataModuleLabel
Definition: SignalToNoise_module.cc:69
lariov::ChannelStatusProvider
Class providing information about the quality of channels.
Definition: ChannelStatusProvider.h:49
geo::GeometryCore::IntersectLines
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.
Definition: GeometryCore.cxx:1320
geo::GeometryCore::View
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Definition: GeometryCore.cxx:601
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
util::absDiff
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
Definition: NumericUtils.h:43
keras_to_tensorflow.int
int
Definition: keras_to_tensorflow.py:69
geo::GeometryCore::PlaneWireToChannel
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
Definition: GeometryCore.cxx:1226
recob::Track::Point_t
tracking::Point_t Point_t
Definition: Track.h:53
geo::GeometryCore::WireEndPoints
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
Definition: GeometryCore.cxx:1243
dune::SignalToNoise::fRawDigitModuleLabel
art::InputTag fRawDigitModuleLabel
Definition: SignalToNoise_module.cc:68
art::fill_ptr_vector
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
cet::sqlite::mean
double mean(sqlite3 *db, std::string const &table_name, std::string const &column_name)
Definition: statistics.cc:16
geo::GeometryCore::WireAngleToVertical
double WireAngleToVertical(geo::View_t view, geo::TPCID const &tpcid) const
Returns the angle of the wires in the specified view from vertical.
Definition: GeometryCore.cxx:843
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
reco_momentum_tuples.x1
x1
Definition: reco_momentum_tuples.py:82
void dune::SignalToNoise::beginJob ( )
overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 436 of file SignalToNoise_module.cc.

437 {
438  art::ServiceHandle<art::TFileService> tfs;
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{
490  art::ServiceHandle<geo::Geometry> geom;
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 }
geo::GeometryCore::AuxDet
AuxDetGeo const & AuxDet(unsigned int const ad=0) const
Returns the specified auxiliary detector.
Definition: GeometryCore.cxx:340
pduneana::line
void line(double t, double *p, double &x, double &y, double &z)
Definition: PDSPAnalyzer_module.cc:4741
geo::GeometryCore::NAuxDets
unsigned int NAuxDets() const
Returns the number of auxiliary detectors.
Definition: GeometryCore.h:4549
event
Event finding and building.
Definition: EventCheater_module.cc:32
SignalToNoise& dune::SignalToNoise::operator= ( SignalToNoise const &  )
delete
SignalToNoise& dune::SignalToNoise::operator= ( SignalToNoise &&  )
delete

Member Data Documentation

double dune::SignalToNoise::cx[93]
private

Definition at line 75 of file SignalToNoise_module.cc.

double dune::SignalToNoise::cy[93]
private

Definition at line 75 of file SignalToNoise_module.cc.

double dune::SignalToNoise::cz[93]
private

Definition at line 75 of file SignalToNoise_module.cc.

TH1D* dune::SignalToNoise::dcos
private

Definition at line 77 of file SignalToNoise_module.cc.

double dune::SignalToNoise::dqdx
private

Definition at line 92 of file SignalToNoise_module.cc.

int dune::SignalToNoise::event
private

Definition at line 88 of file SignalToNoise_module.cc.

art::InputTag dune::SignalToNoise::fCalDataModuleLabel
private

Definition at line 69 of file SignalToNoise_module.cc.

art::InputTag dune::SignalToNoise::fExternalCounterModuleLabel
private

Definition at line 67 of file SignalToNoise_module.cc.

bool dune::SignalToNoise::fIsData
private

Definition at line 71 of file SignalToNoise_module.cc.

trkf::LinFitAlg dune::SignalToNoise::fLinFitAlg
private

Definition at line 73 of file SignalToNoise_module.cc.

art::InputTag dune::SignalToNoise::fRawDigitModuleLabel
private

Definition at line 68 of file SignalToNoise_module.cc.

art::InputTag dune::SignalToNoise::fTrackModuleLabel
private

Definition at line 66 of file SignalToNoise_module.cc.

TTree* dune::SignalToNoise::ftree
private

Definition at line 86 of file SignalToNoise_module.cc.

TH1D* dune::SignalToNoise::hbkg[8][3][4]
private

Definition at line 80 of file SignalToNoise_module.cc.

TH1D* dune::SignalToNoise::hdt
private

Definition at line 84 of file SignalToNoise_module.cc.

TH1D* dune::SignalToNoise::hdx
private

Definition at line 82 of file SignalToNoise_module.cc.

TH2D* dune::SignalToNoise::hphx
private

Definition at line 83 of file SignalToNoise_module.cc.

TH1D* dune::SignalToNoise::hsig[8][3][4]
private

Definition at line 79 of file SignalToNoise_module.cc.

TH1D* dune::SignalToNoise::hstb[8][3][4]
private

Definition at line 81 of file SignalToNoise_module.cc.

TH2D * dune::SignalToNoise::hxz
private

Definition at line 78 of file SignalToNoise_module.cc.

TH2D* dune::SignalToNoise::hyz
private

Definition at line 78 of file SignalToNoise_module.cc.

int dune::SignalToNoise::run
private

Definition at line 87 of file SignalToNoise_module.cc.

int dune::SignalToNoise::tpc
private

Definition at line 89 of file SignalToNoise_module.cc.

int dune::SignalToNoise::trackid
private

Definition at line 91 of file SignalToNoise_module.cc.

int dune::SignalToNoise::wire
private

Definition at line 90 of file SignalToNoise_module.cc.

double dune::SignalToNoise::x
private

Definition at line 93 of file SignalToNoise_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.11