SlicerEva_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file SlicerEva.cxx
3 //
4 // \brief First official DDT analyzer in offline!
5 //
6 // \author ricken@fnal.gov
7 ////////////////////////////////////////////////////////////////////////
15 #include "MCCheater/BackTracker.h"
17 
26 
27 #include "DAQChannelMap/DAQChannelMap.h"
28 #include "RecoBase/Track.h"
29 #include "RecoBase/Cluster.h"
30 #include "Simulation/FLSHit.h"
31 #include "Simulation/FLSHitList.h"
32 
33 
34 #include "TTree.h"
35 
36 
37 #include <string>
38 #include <iostream>
39 #include <cmath>
40 #include <algorithm>
41 #include <vector>
42 #include <fstream>
43 
44 class TTree;
45 
46 namespace ddt {
47  class SlicerEva;
48 }
49 
50 
52 {
53 public:
54  explicit SlicerEva(fhicl::ParameterSet const &pset);
55  virtual ~SlicerEva();
56 
57  void analyze(const art::Event &evt);
58  void beginJob();
59 
60 
61 private:
64  int _ADCthr;
65  int _xMin;
66  int _xMax;
67  int _yMin;
68  int _yMax;
69  int _zMin;
70  int _zMax;
71  int _deltx;
72  int _delty;
73  int _deltz;
74 
75  TTree* _slicetree;
76 
77  float ADC_mean;
78  unsigned long long TDC_min, TDC_max;
79  int X_min, X_max;
80  int Y_min, Y_max;
81  int P_min, P_max;
82  int N_hits, N_sh;
84  int eventID;
85  double beta;
86 };
87 
88 
89 
91  : EDAnalyzer(p),
92  _slicelabel(p.get<std::string>("slice_label")),
93  _sliceinstance(p.get<std::string>("slice_instance")),
94  _ADCthr(p.get<float>("ADC_thr")),
95  _xMin(p.get<float>("x_min")),
96  _xMax(p.get<float>("x_max")),
97  _yMin(p.get<float>("y_min")),
98  _yMax(p.get<float>("y_max")),
99  _deltx(p.get<float>("x_delt")),
100  _delty(p.get<float>("y_delt"))
101 {
102 }
103 
104 //......................................................................
106 
107 //......................................................................
108 
110 {
112  _slicetree = tfs->make<TTree>("slicetree","hit information");
113  _slicetree -> Branch("ADC_mean",&ADC_mean,"ADC_mean/F");
114  _slicetree -> Branch("TDC_min",&TDC_min,"TDC_min/l");
115  _slicetree -> Branch("TDC_max",&TDC_max,"TDC_max/l");
116  _slicetree -> Branch("isMon", &isMon, "isMon/O");
117  _slicetree -> Branch("X_min",&X_min,"X_min/I");
118  _slicetree -> Branch("X_max",&X_max,"X_max/I");
119  _slicetree -> Branch("Y_min",&Y_min,"Y_min/I");
120  _slicetree -> Branch("Y_max",&Y_max,"Y_max/I");
121  _slicetree -> Branch("P_min",&P_min,"P_min/I");
122  _slicetree -> Branch("P_max",&P_max,"P_max/I");
123  _slicetree -> Branch("N_hits",&N_hits,"N_hits/I"); //Number of hits in the slice
124  _slicetree -> Branch("N_sh",&N_sh,"N_sh/I");
125  _slicetree -> Branch("though",&PassThrough,"through/O");
126  _slicetree -> Branch("Event",&eventID,"Event/I");
127  _slicetree -> Branch("beta",&beta,"beta/D");
128  _slicetree -> Branch("contain",&contain,"contain/O");
129  _slicetree -> Branch("wayoff",&wayoff,"wayoff/O");
130  return;
131 }
132 
134 {
135  // contain = false;
136  //Get Truth Info: veto this event if no fast monopole
138  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
139  double t0_MC = 0, t1_MC;
140  bool mon_found = false;
141  std::vector<int> plane_th,cell_th;
142  std::vector<double> time_th;
143 
144  float e_dep = 0;
145 
146  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i){
147  const sim::Particle *p = (*i).second;
148  if (p->PdgCode() == 42) {
149  beta = sqrt((p->P())*(p->P())/((p->P())*(p->P())+ 1.e32) );
150  // if (beta<0.01) return;
151  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());
152  if (flshits.empty()) return;
153  int nFLShits = flshits.size();
154  int starth=0;
155  for (int i =0; i!=nFLShits; ++i) {
156  if (flshits[i].GetPDG()==42) {
157  t0_MC = flshits[i].GetExitT();
158  t1_MC = flshits[i].GetEntryT();
159  starth=i+1;
160  e_dep += flshits[i].GetEdep();
161  plane_th.push_back(flshits[i].GetPlaneID());
162  cell_th.push_back(flshits[i].GetCellID());
163  time_th.push_back(t1_MC);
164 
165  break;
166  }
167 
168  if (starth==nFLShits || starth) return;
169  }
170  for (int h =starth; h!= nFLShits; ++h) {
171  if (flshits[h].GetPDG()!=42) continue;
172  plane_th.push_back(flshits[h].GetPlaneID());
173  cell_th.push_back(flshits[h].GetCellID());
174  time_th.push_back(flshits[h].GetEntryT());
175  e_dep += flshits[h].GetEdep();
176  if (flshits[h].GetExitT()<t0_MC) t0_MC = flshits[h].GetExitT();
177  else if (flshits[h].GetEntryT()>t1_MC) t1_MC = flshits[h].GetEntryT();
178  }
179  }
180  }
181  t0_MC /= 15.625;
182  t1_MC /= 15.625;
183 
184  eventID = e.id().event();
185 
186  //Get the Slice Info:
189  int n_found =0;
190 
191  for (unsigned i = 0; i!= slices->size(); ++i) {
192  isMon = false;
193  novaddt::HitList const &hits = slices->at(i);
194  TDC_min = hits.begin()->TDC().val;
195  TDC_max = TDC_min;
196  P_min = 9999;
197  P_max = -99;
198  X_min = 9999;
199  X_max = -99;
200  Y_min = 9999;
201  Y_max = -99;
202 
203  int ADC_sum = 0;
204  N_hits = hits.size();
205  N_sh = 0;
206 
207  for (auto const hit: hits) {
208  if (TDC_min > hit.TDC().val ) TDC_min = hit.TDC().val;
209  if (TDC_max < hit.TDC().val ) TDC_max = hit.TDC().val;
210  if (P_min > hit.Plane().val) P_min = hit.Plane().val;
211  if (P_max < hit.Plane().val) P_max = hit.Plane().val;
212  if (hit.View().val == daqchannelmap::X_VIEW) {
213  if (X_max < hit.Cell().val) X_max = hit.Cell().val;
214  if (X_min > hit.Cell().val) X_min = hit.Cell().val;
215  }
216  else {
217  if (Y_max < hit.Cell().val) Y_max = hit.Cell().val;
218  if (Y_min > hit.Cell().val) Y_min = hit.Cell().val;
219  }
220 
221  if (hit.ADC().val > _ADCthr) ++N_sh;
222  ADC_sum += hit.ADC().val;
223 
224  if (!mon_found) {
225  for (unsigned hit_th =0 ; hit_th!= plane_th.size(); ++hit_th) {
226  if (hit.Plane().val==plane_th[hit_th] && hit.Cell().val==cell_th[hit_th]) {
227  if (hit.TDC().val< time_th[hit_th]/15.625+500 && hit.TDC().val> time_th[hit_th]/15.625-500)
228  ++n_found;
229  }
230  }
231  if (n_found>0.9*plane_th.size()) {
232  mon_found = true;
233  isMon = true;
234  }
235  }
236  }
237  ADC_mean = ADC_sum/(N_hits+0.);
238  //Decide whether this is contained:
239  if (TDC_min < t0_MC+100 && TDC_max > t1_MC-100) {
240  contain = true;
241  wayoff = false;
242  }
243  else contain = false;
244  if (TDC_min>t1_MC+100) wayoff = true;
245  else wayoff= false;
246  int number =0;
247  if (X_min<_xMin+_deltx) ++number;
248  if (X_max>_xMax-_deltx) ++number;
249  if (Y_min<_yMin+_delty) ++number;
250  if (Y_max>_yMax-_delty) ++number;
251  if (P_min<_zMin+_deltz) ++number;
252  if (P_max>_zMax-_deltz) ++number;
253 
254  if (number>1) PassThrough = true;
255  else PassThrough = false;
256  _slicetree -> Fill();
257  }
258 
259  /*
260  //Here we are implementing some debugging out put file:
261  std::ofstream ofs;
262  ofs.open ("missingMM.dat", std::ofstream::out | std::ofstream::app);
263 
264  if (!mon_found) ofs << plane_th.size() <<" "<<beta<<" "<<e_dep<<std::endl;
265  */
266 
267  return;
268 }
269 
270 
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
list_type::const_iterator const_iterator
DEFINE_ART_MODULE(TestTMapFile)
int TrackId() const
Definition: MCParticle.h:209
unsigned long long TDC_max
unsigned long long TDC_min
void hits()
Definition: readHits.C:15
SlicerEva(fhicl::ParameterSet const &pset)
std::string _sliceinstance
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
double P(const int i=0) const
Definition: MCParticle.h:233
Identifier for the X measuring view of the detector (top)
int evt
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
EventNumber_t event() const
Definition: EventID.h:116
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
void analyze(const art::Event &evt)
std::string _slicelabel
Float_t e
Definition: plot.C:35
EventID id() const
Definition: Event.h:56
enum BeamMode string