All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CRTTrigFilter_module.cc
Go to the documentation of this file.
1 
2 
3 #include "art/Framework/Core/EDFilter.h"
4 #include "art/Framework/Core/ModuleMacros.h"
5 #include "art/Framework/Principal/Event.h"
6 #include "sbnobj/SBND/CRT/CRTData.hh"
7 #include "art_root_io/TFileService.h"
8 #include "TH1.h"
9 #include <bitset>
10 
11 
12 // class CRTTrigFilter : public art::EDFilter(fhicl::ParameterSet const& p) {
13  class CRTTrigFilter : public art::EDFilter {
14  public:
15  explicit CRTTrigFilter(fhicl::ParameterSet const& p);
16  // FlashTimeFilter(FlashTimeFilter const&) = delete;
17  // FlashTimeFilter(FlashTimeFilter&&) = delete;
18  // FlashTimeFilter& operator=(FlashTimeFilter const&) = delete;
19  // FlashTimeFilter& operator=(FlashTimeFilter&&) = delete;
20 
21  //
22  virtual ~CRTTrigFilter() { }
23  virtual bool filter(art::Event& e) override;
24  void reconfigure(fhicl::ParameterSet const& p);
25 
26  private:
27 
28 
29  art::ServiceHandle<art::TFileService> tfs;
30 
31  std::string fCRTStripModuleLabel;
32  std::vector<int> fmodlistUTopL;
33  std::vector<int> fmodlistUBotL;
34  std::vector<int> fmodlistUTopR;
35  std::vector<int> fmodlistUBotR;
36  std::vector<int> fmodlistDTopL;
37  std::vector<int> fmodlistDBotL;
38  std::vector<int> fmodlistDTopR;
39  std::vector<int> fmodlistDBotR;
40  bool fstripmatch;
41  int fstripshift;
42  int fedgecut;
43  float fTimeCoinc;
44  float fADCthresh;
45 
46 
47  TH1F *hits;
48  TH1F *trigt;
49  TH1F *trig;
50 
51  };
52 
53  CRTTrigFilter::CRTTrigFilter(fhicl::ParameterSet const& p): EDFilter{p}
54  {
55 
56 
57  // Create histograms
58  hits = tfs->make<TH1F>("hits","hits",14,39.5,63.5);
59  // hits->GetXaxis()->SetTitle("Track Length (cm)");
60  trigt = tfs->make<TH1F>("trigt","trigt",200,-3000,3000); // trig time in us
61  trig = tfs->make<TH1F>("trig","trig",2,-0.5,1.5);
62 
63  this->reconfigure(p);
64  }
65 
66  void CRTTrigFilter::reconfigure(fhicl::ParameterSet const& p)
67  {
68 
69  fCRTStripModuleLabel = p.get< std::string >("CRTStripModuleLabel","crt");
70  fmodlistUTopL = p.get<std::vector<int>>("ModuleListUpstreamTopLeft");
71  fmodlistUBotL = p.get<std::vector<int>>("ModuleListUpstreamBotLeft");
72  fmodlistUTopR = p.get<std::vector<int>>("ModuleListUpstreamTopRight");
73  fmodlistUBotR = p.get<std::vector<int>>("ModuleListUpstreamBotRight");
74  fmodlistDTopL = p.get<std::vector<int>>("ModuleListDownstreamTopLeft");
75  fmodlistDBotL = p.get<std::vector<int>>("ModuleListDownstreamBotLeft");
76  fmodlistDTopR = p.get<std::vector<int>>("ModuleListDownstreamTopRight");
77  fmodlistDBotR = p.get<std::vector<int>>("ModuleListDownstreamBotRight");
78  fstripmatch = p.get<bool>("RequireStripMatch",true);
79  fstripshift = p.get<int>("AllowedStripShift",4);
80  fedgecut = p.get<int>("EdgeCutInStrips",12);
81  if (fedgecut>31) fedgecut=31;
82  fTimeCoinc = p.get<float>("StripTimeCoincidence",0.2);
83  fADCthresh = p.get<float>("ADCthresh",500.0);
84  }
85 
86  bool CRTTrigFilter::filter(art::Event& e)
87  {
88 
89  bool KeepMe = false;
90  int event = e.id().event();
91  if (event%1000==0) std::cout << "event " << event << std::endl;
92  // unsigned long long planeUtop, planeUbot, planeDtop, planeDbot;
93  uint64_t planeUtopL, planeUbotL, planeDtopL, planeDbotL;
94  uint64_t planeUtopR, planeUbotR, planeDtopR, planeDbotR;
95  uint64_t planeUpStL, planeDownStL;
96  uint64_t planeUpStR, planeDownStR;
97 
98  // uint64 planeUpSt,planeDownSt;
99 
100  int nstr=0;
101  art::Handle<std::vector<sbnd::crt::CRTData> > crtStripListHandle;
102  std::vector<art::Ptr<sbnd::crt::CRTData> > striplist;
103  if (e.getByLabel(fCRTStripModuleLabel, crtStripListHandle)) {
104  // if (e.getByLabel("crt", crtStripListHandle)) {
105  art::fill_ptr_vector(striplist, crtStripListHandle);
106  nstr = striplist.size();
107  }
108  // std::cout << "number of crt strips " << nstr << std::endl;
109 
110  bool trigKeep = false;
111  uint64_t edgecutL,edgecutR;
112  edgecutR = ~((uint64_t)(pow(2,fedgecut)-1) | ((uint64_t)(pow(2,32)-1)<< 32));
113  edgecutL = ~(((uint64_t)(pow(2,fedgecut)-1) << (32-fedgecut))| ((uint64_t)(pow(2,32)-1) <<32));
114  // std::cout << "left " << std::bitset<64>(edgecutL) << " right " << std::bitset<64>(edgecutR) << std::endl;
115  // left 0000000000000000000000000000000000000000000011111111111111111111
116  // right 0000000000000000000000000000000011111111111111111111000000000000
117 
118  float stripwidth = 11.2; //cm
119  float driftvel = 0.16; //cm/us
120 
121  // set cut limits on drift window time expectation of track from CRT hit.
122  float rwcut = fedgecut*stripwidth/driftvel;
123  float rwcutlow = -200.0 - rwcut;
124  float rwcuthigh = 1500.00 + rwcut;
125 
126  for (int i = 0; i<nstr-2; i+=2){
127  if ((striplist[i]->ADC()+striplist[i+1]->ADC())>fADCthresh) {
128  uint32_t chan1 = striplist[i]->Channel();
129  int strip1 = (chan1 >> 1) & 15;
130  int module1 = (chan1>> 5);
131  uint32_t ttime1 = striplist[i]->T0();
132  float ctime1 = ttime1/16.;
133  if (ttime1 > 2147483648) {
134  ctime1 = ((ttime1-4294967296)/16.);
135  }
136 
137  //
138  for (int j = i+2; j<nstr; j+=2){
139  if ((striplist[j]->ADC()+striplist[j+1]->ADC())>fADCthresh) {
140  bool match = false;
141  uint32_t chan2 = striplist[j]->Channel();
142  int strip2 = (chan2 >> 1) & 15;
143  int module2 = (chan2>> 5);
144  //
145  uint32_t ttime2 = striplist[j]->T0();
146  // T0 in units of ticks, but clock frequency (16 ticks = 1 us) is wrong.
147  // ints were stored as uints, need to patch this up
148  float ctime2 = ttime2/16.;
149  if (ttime2 > 2147483648) {
150  ctime2 = ((ttime2-4294967296)/16.);
151  }
152  //
153  float diff = abs(ctime1-ctime2);
154  // if (diff<=fTimeCoinc && ctime1>-1300 && ctime1<1300) {
155  if (diff<=fTimeCoinc ) {
156 
157 
158  planeUtopL=0; planeUbotL=0; planeDtopL=0; planeDbotL=0;
159  planeUtopR=0; planeUbotR=0; planeDtopR=0; planeDbotR=0;
160  for (size_t im=0;im<fmodlistUTopL.size();++im) {
161  size_t iswap = fmodlistUTopL.size()-im-1;
162  if (module1==fmodlistUTopL[im]) planeUtopL+=(1 << ((15-strip1)+iswap*16));
163  if (module2==fmodlistUTopL[im]) planeUtopL+=(1 << ((15-strip2)+iswap*16));
164  }
165  for (size_t im=0;im<fmodlistUBotL.size();++im) {
166  size_t iswap = fmodlistUBotL.size()-im-1;
167  if (module1==fmodlistUBotL[im]) planeUbotL+=(1 << ((15-strip1)+iswap*16));
168  if (module2==fmodlistUBotL[im]) planeUbotL+=(1 << ((15-strip2)+iswap*16));
169  }
170  for (size_t im=0;im<fmodlistUTopR.size();++im) {
171  size_t iswap = fmodlistUTopR.size()-im-1;
172  if (module1==fmodlistUTopR[im]) planeUtopR+=(1 << ((15-strip1)+iswap*16));
173  if (module2==fmodlistUTopR[im]) planeUtopR+=(1 << ((15-strip2)+iswap*16));
174  }
175  for (size_t im=0;im<fmodlistUBotR.size();++im) {
176  size_t iswap = fmodlistUBotR.size()-im-1;
177  if (module1==fmodlistUBotR[im]) planeUbotR+=(1 << ((15-strip1)+iswap*16));
178  if (module2==fmodlistUBotR[im]) planeUbotR+=(1 << ((15-strip2)+iswap*16));
179  }
180  for (size_t im=0;im<fmodlistDTopL.size();++im) {
181  size_t iswap = fmodlistDTopL.size()-im-1;
182  if (module1==fmodlistDTopL[im]) planeDtopL+=(1 << ((15-strip1)+iswap*16));
183  if (module2==fmodlistDTopL[im]) planeDtopL+=(1 << ((15-strip2)+iswap*16));
184  }
185  for (size_t im=0;im<fmodlistDBotL.size();++im) {
186  size_t iswap = fmodlistDBotL.size()-im-1;
187  if (module1==fmodlistDBotL[im]) planeDbotL+=(1 << ((15-strip1)+iswap*16));
188  if (module2==fmodlistDBotL[im]) planeDbotL+=(1 << ((15-strip2)+iswap*16));
189  }
190  for (size_t im=0;im<fmodlistDTopR.size();++im) {
191  size_t iswap = fmodlistDTopR.size()-im-1;
192  if (module1==fmodlistDTopR[im]) planeDtopR+=(1 << ((15-strip1)+iswap*16));
193  if (module2==fmodlistDTopR[im]) planeDtopR+=(1 << ((15-strip2)+iswap*16));
194  }
195  for (size_t im=0;im<fmodlistDBotR.size();++im) {
196  size_t iswap = fmodlistDBotR.size()-im-1;
197  if (module1==fmodlistDBotR[im]) planeDbotR+=(1 << ((15-strip1)+iswap*16));
198  if (module2==fmodlistDBotR[im]) planeDbotR+=(1 << ((15-strip2)+iswap*16));
199  }
200 
201  planeUpStL = planeUtopL | planeUbotL;
202  planeUpStR = planeUtopR | planeUbotR;
203  planeDownStL = planeDtopL | planeDbotL;
204  planeDownStR = planeDtopR | planeDbotR;
205  if ((planeUpStR & edgecutR) && (planeDownStR & edgecutR)) {
206 
207  if ( planeUpStR & planeDownStR ) match=true;
208  for (int is=1;is<=fstripshift;++is) {
209  uint64_t temp; temp=(planeUpStR << is);
210  if (temp & planeDownStR) match=true;
211  temp=(planeDownStR << is);
212  if (temp & planeUpStR) match=true;
213  }
214  }
215  if ((planeUpStL & edgecutL) && (planeDownStL & edgecutL)) {
216  if ( planeUpStL & planeDownStL ) match=true;
217  for (int is=1;is<=fstripshift;++is) {
218  uint64_t temp; temp=(planeUpStL << is);
219  if (temp & planeDownStL) match=true;
220  temp=(planeDownStL << is);
221  if (temp & planeUpStL) match=true;
222  }
223  }
224  if (match) {
225  // std::cout << " before drift window cut " << std::endl;
226  // std::cout << " module and strip 1 " << module1 << " " << strip1 << std::endl;
227  // std::cout << " module and strip 2 " << module2 << " " << strip2 << std::endl;
228  // std::cout << " upstream " << std::bitset<64>(planeUpStR) <<
229  // std::endl << "downstream " << std::bitset<64>(planeDownStR) <<
230  // std::endl;
231  // std::cout << " ----------------------------------" << std::endl;
232  float xpos; // in cm
233  if (module1<module2) xpos=((int(module1/2)-20)*16+strip1+0.5)*stripwidth-358.4;
234  else xpos=((int(module2/2)-20)*16+strip2+0.5)*stripwidth-358.4;
235  float dtime; // in us
236  if (xpos>0) dtime = ctime1+((200.0-xpos)/driftvel);
237  else dtime = ctime1+((200.0+xpos)/driftvel);
238  if (dtime<rwcutlow || dtime>rwcuthigh ) match=false;
239  // std::cout << "Event " << event << " " << module1 << " " << xpos << " " << ctime1 << " " << dtime << std::endl;
240  }
241  }
242  if (match) {
243  // std::cout << "Found One! Event " << event <<
244  // " m/s1 m/s2 " << module1 << " " << strip1 << " " <<
245  // module2 << " " << strip2 << " " << ctime1 << " " <<ctime2 << std::endl;
246  KeepMe=true;
247  }
248  }} // end loop over second strip
249  }} // end loop over first strip
250  if (trigKeep) trig->Fill(1.0); else trig->Fill(0.0);
251  return KeepMe;
252 
253  }
254 
255  // A macro required for a JobControl module.
256  DEFINE_ART_MODULE(CRTTrigFilter)
257 
258 
CRTTrigFilter::tfs
art::ServiceHandle< art::TFileService > tfs
Definition: CRTTrigFilter_module.cc:29
CRTTrigFilter::reconfigure
void reconfigure(fhicl::ParameterSet const &p)
Definition: CRTTrigFilter_module.cc:66
p
pdgs p
Definition: selectors.fcl:22
CRTTrigFilter::fmodlistUBotL
std::vector< int > fmodlistUBotL
Definition: CRTTrigFilter_module.cc:33
reconfigure
for pfile in ack l reconfigure(.*) override"` do echo "checking $
Definition: no_reconfigure_override.sh:3
CRTTrigFilter::fmodlistUBotR
std::vector< int > fmodlistUBotR
Definition: CRTTrigFilter_module.cc:35
CRTTrigFilter::fmodlistDBotL
std::vector< int > fmodlistDBotL
Definition: CRTTrigFilter_module.cc:37
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
CRTTrigFilter::filter
virtual bool filter(art::Event &e) override
Definition: CRTTrigFilter_module.cc:86
CRTTrigFilter::fmodlistDTopL
std::vector< int > fmodlistDTopL
Definition: CRTTrigFilter_module.cc:36
CRTTrigFilter::fCRTStripModuleLabel
std::string fCRTStripModuleLabel
Definition: CRTTrigFilter_module.cc:31
CRTTrigFilter::fmodlistDTopR
std::vector< int > fmodlistDTopR
Definition: CRTTrigFilter_module.cc:38
CRTTrigFilter::fmodlistUTopL
std::vector< int > fmodlistUTopL
Definition: CRTTrigFilter_module.cc:32
CRTTrigFilter::CRTTrigFilter
CRTTrigFilter(fhicl::ParameterSet const &p)
Definition: CRTTrigFilter_module.cc:53
CRTTrigFilter::fmodlistUTopR
std::vector< int > fmodlistUTopR
Definition: CRTTrigFilter_module.cc:34
e
do i e
Definition: updateVersionList.sh:68
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
CRTTrigFilter::fmodlistDBotR
std::vector< int > fmodlistDBotR
Definition: CRTTrigFilter_module.cc:39
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
Generated by   doxygen 1.8.5