ToFPositionRecoAnalysis_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 /// \file ToFPositionRecoAnalysis_module.cc
3 /// \module analyzer
4 /// \brief Experiment with position reconstruction
5 /// \author John Rabaey, University of Dallas, <jrabaey@udallas.edu>
6 /// \date March 2019
7 ////////////////////////////////////////////////////////////////////////////
8 
9 // framework
13 #include "fhiclcpp/ParameterSet.h"
22 
23 // nova
24 #include "RawData/RawBeamline.h"
28 #include "BeamlineRecoBase/ToF.h"
32 
33 // ROOT
34 #include "TH1D.h"
35 #include "TH1I.h"
36 #include "TH2D.h"
37 #include "TCanvas.h"
38 #include "TGraph.h"
39 
40 // stl
41 #include <iostream>
42 #include <algorithm>
43 #include <fstream>
44 
45 //a whole bunch of histograms?
46 bool dumpPosns = false;
47 
48 //main histogram window constants
49 double MIN_TOF = 40;
50 double MAX_TOF = 47;
51 double N_BINS = 100;
52 
53 double MIN_DT = -5;
54 double MAX_DT = 5;
55 
56 double MAX_SPREAD = 4;
57 double SPREAD_CUTOFF = 4;
58 
59 bool POSN_ZOOM_IN = true;
60 
61 //position reconstruction constants
62 bool VARY_SIDE = false;
63 double SCINT_SIDE_LENGTH = 2.5; // In ns // use if vary_side is false
64 double MIN_SCINT_SIDE = 0; //use if vary_side is true
65 double MAX_SCINT_SIDE = 4.0;
66 double SCINT_SIDE_STEP = .1;
67 
68 
69 //position reconstruction dump binning
70 int POSN_N_BINS = 40;
71 
72 enum ToFCounter {
77 };
78 
79 //channel numbers
80 const int nChannels = 4;
81 std::map<ToFCounter,std::vector<unsigned int>> channels = {
82  { kToFCounterDS, {15,14,13,12} },
83  { kToFCounterUS, {8, 11,10,9 } },
84  { kToFCounterDSSiPM,{7, 6, 5, 4 } }
85 };
86 std::vector<unsigned int> requiredChannels = {};// channels[kToFCounterDSSiPM];
87 
88 std::vector<ToFCounter> activeToFCounters = {
92 };
93 
94 std::map<ToFCounter,TString> tofCounterStrings = {
95  {kToFCounterUS, "USToF__"},
96  {kToFCounterDS, "DSToF__"},
97  {kToFCounterDSSiPM, "SiPMToF"},
98  {kToFCounterError, "ErrorTF"}
99 };
100 
101 std::map<ToFCounter,double> tofCounterDelays = {
102  { kToFCounterUS, 21.0 },
103  { kToFCounterDS, 0 },
104  { kToFCounterDSSiPM, 0 }
105 };
106 
107 // -----------------------------------------------------------------------
108 namespace beamlinereco {
110 
111  public:
112 
114 
115  void reconfigure(const fhicl::ParameterSet& pset);
116  void beginJob() override;
117  void analyze(const art::Event& evt);
118  void endJob() override;
119 
120  private:
121 
123 
124  std::unordered_map<unsigned int, std::string> channelMap;
125 
128 
129  std::vector<std::pair<unsigned int, std::string>> fChannels;
131 
132  TH1D* hToF;
133  std::map<int, std::map<int, TH1D*>> hdtByChannel;
134  std::map<ToFCounter,TH1D*> hFirstHits;
135  std::map<ToFCounter,std::vector<TH2D*>> hPosnRecoDump;
136  std::map<ToFCounter,std::vector<TH2D*>> hPosn;
137  std::map<ToFCounter,std::vector<TH1D*>> hPosnSpread;
138 
139  std::map<ToFCounter,std::vector<int>> nPhysicalPosns;
140  std::map<ToFCounter,std::vector<int>> nNonphysicalPosns;
141  std::map<ToFCounter,int> nFourHitEvents;
142  std::map<ToFCounter,int> nEvents;
143  std::map<ToFCounter,int> nHighSpreadPosns;
144 
146 
147  // services
149  };
150 
152 
153 }
154 
155 // -----------------------------------------------------------------------
157  this->reconfigure(pset);
158 }
159 
160 // -----------------------------------------------------------------------
162  fRecoToFLabel = pset.get<art::InputTag>("RecoToFLabel");
163  fToFDigitLabel = pset.get<art::InputTag>("ToFDigitLabel");
164 
165 
166  fChannels = pset.get<std::vector<std::pair<unsigned int, std::string>>>("Channels");
167  fSamplingInterval = pset.get<double>("TimeSamplingInterval");
168 
169  channelMap.clear();
170  for (auto ch : fChannels) { channelMap[ch.first] = ch.second; }
171 }
172 
173 // -----------------------------------------------------------------------
175 
176  //set up directories, main ToF distributions
177  this->hToF = tfs->make<TH1D>("hToF", "Time of Flight", N_BINS, MIN_TOF, MAX_TOF);
178 
179  //Do everything three times, once for each ToF
180  for(auto tofc : activeToFCounters) {
181  art::TFileDirectory tofDir = tfs->mkdir(tofCounterStrings[tofc].Data());
182 
183  //set up dt distributions
184  art::TFileDirectory dtDir = tofDir.mkdir("dt_distributions");
185 
186  for(int i = 0; i < nChannels; i++) {
187  for(int j = i+1; j < nChannels; j++) {
188  hdtByChannel[channels[tofc][i]][channels[tofc][j]] =
189  dtDir.make<TH1D>(Form("h_dt_ch%d_%d", channels[tofc][i], channels[tofc][j]),
190  Form("dt for channels %d and %d", channels[tofc][i], channels[tofc][j]),
191  N_BINS, MIN_DT, MAX_DT);
192  }//for j
193  }//for i
194 
195  //set up plot for the first channel in each hit group
196  hFirstHits[tofc] = tofDir.make<TH1D>("first_hits", "Number of times each channel is first in a cluster",4,0,4);
197 
198  //set up position distribution
199  if(VARY_SIDE) {
200  //better put things in a folder
201  art::TFileDirectory posnstuff = tofDir.mkdir("hPosn");
202  //vector needs to contain one double per side length
203  for(double d = MIN_SCINT_SIDE; d < MAX_SCINT_SIDE; d += SCINT_SIDE_STEP) {
204  this->hPosn[tofc].push_back(posnstuff.make<TH2D>(Form("hPosnS%.3f",d),
205  Form("Hit position with side length constant %.3f", d),
206  N_BINS,
207  (POSN_ZOOM_IN ? 0 : -d),
208  (POSN_ZOOM_IN ? d : 2*d),
209  N_BINS,
210  (POSN_ZOOM_IN ? 0 : -d),
211  (POSN_ZOOM_IN ? d : 2*d)));
212  this->hPosnSpread[tofc].push_back(posnstuff.make<TH1D>(Form("hPosnSpreadS%.3f", d),
213  Form("Position spread with side length constant %.3f", d),
214  N_BINS, 0, MAX_SPREAD));
215  nPhysicalPosns[tofc].push_back(0);
216  nNonphysicalPosns[tofc].push_back(0);
217  } //for d
218  } //if
219 
220  else {//(!VARY_SIDE)
221  //vector only needs one thing
222  this->hPosn[tofc].push_back(tofDir.make<TH2D>("hPosn", "Position",
223  N_BINS,
226  N_BINS,
227  (POSN_ZOOM_IN ? 0 : -SCINT_SIDE_LENGTH),
228  (POSN_ZOOM_IN ? SCINT_SIDE_LENGTH : 2*SCINT_SIDE_LENGTH)));
229  this->hPosnSpread[tofc].push_back(tofDir.make<TH1D>("hPosnSpread", "Position spread",
230  N_BINS, 0, MAX_SPREAD));
231  nPhysicalPosns[tofc].push_back(0);
232  nNonphysicalPosns[tofc].push_back(0);
233  } //else
234 
235  //set draw options for 2d histograms
236  for(TH2D* posnHist : hPosn[tofc]) {
237  posnHist->SetOption("colz");
238  }
239 
240  //set the useful counter
241  nEvents[tofc] = 0;
242  nHighSpreadPosns[tofc] = 0;
243  nFourHitEvents[tofc] = 0;
244 
245  }//for tofc
246 
248 }
249 
250 // -----------------------------------------------------------------------
252 
253  // Get ToFs for each reco version
255  std::vector<art::Ptr<brb::ToF>> tofs;
256  if (e.getByLabel(fRecoToFLabel, tofHandle)) {
257  art::fill_ptr_vector(tofs, tofHandle);
258  }
259  unsigned int nToFs = tofs.size();
260 
261  // Get time of flight for each ToF
262  std::vector<double> flightTimes;
263  for(unsigned int i = 0; i < nToFs; i++) {
264  flightTimes.push_back(tofs[i]->Time());
265  }
266 
267  //fill tof histogram
268  for(unsigned int i = 0; i < nToFs; i++) {
269  hToF->Fill(flightTimes[i]);
270  }
271 
272 
273  // Get associated digits for each reco version
275  std::vector<art::Ptr<brb::BeamlineDigit>> tofDigits;
276  if (e.getByLabel(fToFDigitLabel, tofDigitHandle)) {
277  art::fill_ptr_vector(tofDigits, tofDigitHandle);
278  }
279 
280  // Retrieve start times by channel
281  std::unordered_map<int,double> startTimeByChannel;
282  //first, initialize all required channels to -1. This flag will allow us to skip the event
283  // if the channel doesn't get a good start time value
284  for(auto tofc : activeToFCounters) {
285  for(int ch : channels[tofc]) {
286  startTimeByChannel[ch] = -1;
287  }//for ch
288  }//for tofc
289 
290  for(auto d : tofDigits){
291  //get start time by channel
292  startTimeByChannel[fChannelMap->OnlineDigitChannel(d->ChannelID())] = d->StartTimeInNanoSec();
293  }//for
294 
295  //skip the event if it doesn't have all required channels
296  for(int ch : requiredChannels) {
297  if(startTimeByChannel[ch] <= 0) {
298  return;
299  }//if
300  }//for
301 
302  //debugging
303  //std::cout << "Got start times for run " << e.run() << " subrun " << e.subRun() << " event " << e.id().event() << "\n";
304 
305  //everything else needs to be wrapped in a big for loop and executed once for each tof counter
306  for(auto tofc : activeToFCounters) {
307  //*****************THE DT DISTRIBUTIONS*******************
308  for (int i=0; i<4; i++) {
309  for (int j=i+1; j<4; j++) {
310  if(startTimeByChannel[channels[tofc][i]] > 0 && startTimeByChannel[channels[tofc][j]] > 0) {
311  double dt = startTimeByChannel[channels[tofc][i]] - startTimeByChannel[channels[tofc][j]];
312  hdtByChannel[channels[tofc][i]][channels[tofc][j]]->Fill(dt);
313  }//if
314  }//for i
315  }//for j
316 
317  //debugging
318  //std::cout << "Made dt distributions for run " << e.run() << " subrun " << e.subRun() << " event " << e.id().event() << "\n";
319 
320  //************THE FIRST CHANNEL HIT DISTRIBUTIONS*************
321  bool fourHits = true;
322  for(int i=0; i<4; i++) {
323  fourHits = fourHits && (startTimeByChannel[channels[tofc][i]] > 0);
324  }
325  int firstBin = 0;
326  for(int i=0; i<4; i++) {
327  firstBin = startTimeByChannel[channels[tofc][i]] < startTimeByChannel[channels[tofc][firstBin]] ? i : firstBin;
328  } //for i
329  if(fourHits) {
330  hFirstHits[tofc]->Fill(firstBin);
331  }
332 
333  //debugging
334  //std::cout << "Made first channel distributions for run " << e.run() << " subrun " << e.subRun() << " event " << e.id().event() << "\n";
335 
336  //******************POSITION RECONSTRUCTION DUMP**************
337  //only if hits recorded on all four corners
338  if(fourHits) {
339  if(dumpPosns) {
340  //make histograms
341  TH2D* eventPosnDump
342  = new TH2D(posnfinder->GetSummaryHist(Form("run%d_subrun%d_event%d_posndump", e.run(), e.subRun(), e.id().event()),
343  Form("run%d_subrun%d_event%d_posndump", e.run(), e.subRun(), e.id().event()),
344  startTimeByChannel[channels[tofc][0]],
345  startTimeByChannel[channels[tofc][1]],
346  startTimeByChannel[channels[tofc][2]],
347  startTimeByChannel[channels[tofc][3]],
349  hPosnRecoDump[tofc].push_back(eventPosnDump);
350  }
351 
352  if(VARY_SIDE) {
353  int i = 0; //got to keep track in two places...
354  for(double s = MIN_SCINT_SIDE; s < MAX_SCINT_SIDE; s += SCINT_SIDE_STEP) {
355  double spread = posnfinder->Get4CornerSpread(startTimeByChannel[channels[tofc][0]],
356  startTimeByChannel[channels[tofc][1]],
357  startTimeByChannel[channels[tofc][2]],
358  startTimeByChannel[channels[tofc][3]],
359  s);
360  hPosnSpread[tofc][i]->Fill(spread);
361 
362  posn_t loc = posnfinder->Get4CornerPosition(startTimeByChannel[channels[tofc][0]],
363  startTimeByChannel[channels[tofc][1]],
364  startTimeByChannel[channels[tofc][2]],
365  startTimeByChannel[channels[tofc][3]],
366  s);
367  if(spread < SPREAD_CUTOFF) {
368  hPosn[tofc][i]->Fill(loc.x,loc.y);
369  bool phys = posnfinder->IsPhysical(loc,s);
370  nPhysicalPosns[tofc][i] += phys ? 1 : 0;
371  nNonphysicalPosns[tofc][i] += phys ? 0 : 1;
372  }//if
373  i++;
374  }//for
375  }//if
376 
377  else { //(!VARY_SIDE)
378  double t1 = startTimeByChannel[channels[tofc][0]];
379  double t2 = startTimeByChannel[channels[tofc][1]];
380  double t3 = startTimeByChannel[channels[tofc][2]];
381  double t4 = startTimeByChannel[channels[tofc][3]];
382  double spread = posnfinder->Get4CornerSpread(t1,t2,t3,t4,SCINT_SIDE_LENGTH);
383  hPosnSpread[tofc][0]->Fill(spread);
384 
386  if(spread < SPREAD_CUTOFF) {
387  hPosn[tofc][0]->Fill(loc.x,loc.y);
388 
389  bool phys = posnfinder->IsPhysical(loc,SCINT_SIDE_LENGTH);
390  nPhysicalPosns[tofc][0] += phys ? 1 : 0;
391  nNonphysicalPosns[tofc][0] += phys ? 0 : 1;
392 
393  }//if
394  else {
395  std::cout << Form("High spread for run %d subrun %d event %d with hits at %.2f,%.2f,%.2f,%.2f and reco posn %.2f,%.2f: %.2f\n", e.run(), e.subRun(), e.id().event(), t1, t2, t3, t4, loc.x, loc.y, spread);
396  nHighSpreadPosns[tofc]++;
397  }//else
398  }//else
399 
400  nFourHitEvents[tofc]++;
401  }//if fourhits
402 
403 
404  //keep count
405  nEvents[tofc]++;
406 
407 
408  //debugging
409  //std::cout << "Made position histograms for run " << e.run() << " subrun " << e.subRun() << " event " << e.id().event() << "\n";
410  }//for tofc
411 }//analyze(event)
412 
414  if(dumpPosns) {
415  art::TFileDirectory posnDumpDir = tfs->mkdir("Position_reco_dump");
416  for(auto tofc : activeToFCounters) {
417  //send the posn reco dump to a folder
418  TString name = "Posns_" + tofCounterStrings[tofc];
419  art::TFileDirectory posnRecoDumpDir = posnDumpDir.mkdir(name.Data());
420  for(auto eventHist : hPosnRecoDump[tofc]) {
421  posnRecoDumpDir.make<TH2D>(*eventHist);
422  }//for eventHist
423  }//for tofc
424  }//if dumpPosns
425 
426  //prep file
427  std::ofstream dtInfo;
428  dtInfo.open("dtInfo.txt");
429  for(auto tofc : activeToFCounters) {
430  // Write all ToF dt information to a text file for easy reading and viewing
431  for (int i=0; i<4; i++) {
432  for (int j=i+1; j<4; j++) {
433  //calculate values
434  double mean = hdtByChannel[channels[tofc][i]][channels[tofc][j]]->GetMean(1); //mean of the histogram along the x=1 axis
435  double rms = hdtByChannel[channels[tofc][i]][channels[tofc][j]]->GetRMS(1); //standard deviation
436 
437  //write to file
438  dtInfo << std::string(40,'-') << "\n";
439  dtInfo << Form("| %2d | %2d | %7.4f | %6.4f |\n", channels[tofc][i], channels[tofc][j], mean, rms);
440  } //for i
441  } //for j
442  } //for tofc
443  dtInfo.close();
444 
445  for(auto tofc : activeToFCounters) {
446  //print the number of physical and all posns
447  if(VARY_SIDE) {
448  int i = 0;
449  for(double s = MIN_SCINT_SIDE; s < MAX_SCINT_SIDE; s += SCINT_SIDE_STEP) {
450  std::cout << Form("In " + tofCounterStrings[tofc] + " There were %d physical posns ", nPhysicalPosns[tofc][i])
451  << Form("for side = %.3f\n", s);
452  std::cout << Form("In " + tofCounterStrings[tofc] + " There were %d nonphysical posns ", nNonphysicalPosns[tofc][i])
453  << Form("for side = %.3f\n", s);
454  i++;
455  }//for
456  }//if
457  else { //!VARY_SIDE
458  std::cout << Form("In " + tofCounterStrings[tofc] + " There were %d physical posns.\n", nPhysicalPosns[tofc][0]);
459  std::cout << Form("In " + tofCounterStrings[tofc] + " There were %d nonphysical posns.\n", nNonphysicalPosns[tofc][0]);
460  }//else
461 
462  std::cout << Form("In total, for " + tofCounterStrings[tofc] + " there were %d events with spread above %.3f.\n", nHighSpreadPosns[tofc], SPREAD_CUTOFF);
463  std::cout << Form("In total, for " + tofCounterStrings[tofc] + " there were %d events with four hits.\n", nFourHitEvents[tofc]);
464  std::cout << Form("In total, for " + tofCounterStrings[tofc] + " there were %d events.\n", nEvents[tofc]);
465 
466  //redo the labels for the first hits histograms
467  //and set a nice range
468  for(int i = 0; i < 4; i++) {
469  hFirstHits[tofc]->GetXaxis()->SetBinLabel(i+1, Form("ch%d",channels[tofc][i]));
470  hFirstHits[tofc]->GetYaxis()->SetRangeUser(0, 1.2*hFirstHits[tofc]->GetBinContent(hFirstHits[tofc]->GetMaximumBin()));
471  hFirstHits[tofc]->SetOption("bar1");
472  }//for i
473 
474  }//for tofc
475 }//endJob
const XML_Char * name
Definition: expat.h:151
SubRunNumber_t subRun() const
Definition: Event.h:72
std::map< ToFCounter, double > tofCounterDelays
TH2D GetSummaryHist(TString name, TString title, double tA, double tB, double tC, double tD, double s)
std::map< ToFCounter, std::vector< TH2D * > > hPosn
double SCINT_SIDE_LENGTH
art::ServiceHandle< beamlineutil::BeamlineChannelMap > fChannelMap
std::map< ToFCounter, std::vector< TH2D * > > hPosnRecoDump
std::map< ToFCounter, TString > tofCounterStrings
TFileDirectory mkdir(std::string const &dir, std::string const &descr="")
std::map< ToFCounter, std::vector< int > > nPhysicalPosns
const int nChannels
DEFINE_ART_MODULE(TestTMapFile)
unsigned int OnlineDigitChannel(ChannelID channel) const
Online digiziter channel number for this offline ChannelID.
std::vector< unsigned int > requiredChannels
const XML_Char * s
Definition: expat.h:262
Definition: Cand.cxx:23
std::map< ToFCounter, std::vector< unsigned int > > channels
double Get4CornerSpread(double tA, double tB, double tC, double tD, double s)
Encapsulation of reconstructed digitizer &#39;hits&#39;. Used for ToF PMTs and SiPMs, and Cherenkov and Muon ...
std::vector< std::pair< unsigned int, std::string > > fChannels
T get(std::string const &key) const
Definition: ParameterSet.h:231
Encapsulation of reconstructed Time-of-Flight (ToF) information. Part of beamline reconstruction for ...
int evt
bool IsPhysical(posn_t p, double s)
Float_t d
Definition: plot.C:236
Encapsulation of reconstructed Wire Chamber track. Part of beamline reconstruction for NOvA test beam...
const double j
Definition: BetheBloch.cxx:29
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double t2
OStream cout
Definition: OStream.cxx:6
std::map< int, std::map< int, TH1D * > > hdtByChannel
posn_t Get4CornerPosition(double tA, double tB, double tC, double tD, double s)
EventNumber_t event() const
Definition: EventID.h:116
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
art::ServiceHandle< art::TFileService > tfs
Raw data definitions for beamline data used in NOvA test beam experiment.
std::map< ToFCounter, std::vector< int > > nNonphysicalPosns
ToFPositionRecoAnalysis(const fhicl::ParameterSet &pset)
void reconfigure(const fhicl::ParameterSet &pset)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Encapsulation of reconstructed PID information from detectors in the beamline (ToF, WCs, Cherenkov). Part of beamline reconstruction for NOvA test beam.
double SCINT_SIDE_STEP
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
std::vector< ToFCounter > activeToFCounters
EventID id() const
Definition: Event.h:56
Channel mapping service which may be used to interpret the channels which are read out by the various...
std::unordered_map< unsigned int, std::string > channelMap
Encapsulation of reconstructed track in the muon stack downstream of test beam detector. Part of beamline reconstruction for NOvA test beam.
std::map< ToFCounter, std::vector< TH1D * > > hPosnSpread