All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.cc
Go to the documentation of this file.
1 #include "CRTHitRecoAlg.h"
2 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
3 #include <algorithm>
4 using namespace icarus::crt;
5 
6 //----------------------------------------------------------------------
7 CRTHitRecoAlg::CRTHitRecoAlg(const Config& config){
8  this->reconfigure(config);
9  fChannelMap = art::ServiceHandle<icarusDB::IICARUSChannelMap const>{}.get();
10  fGeometryService = lar::providerFrom<geo::Geometry>();
11  fCrtutils = new CRTCommonUtils();
12 }
13 
14 //---------------------------------------------------------------------
15 CRTHitRecoAlg::CRTHitRecoAlg(){
16  fChannelMap = art::ServiceHandle<icarusDB::IICARUSChannelMap const>{}.get();
17  fGeometryService = lar::providerFrom<geo::Geometry>();
18  fCrtutils = new CRTCommonUtils();
19 }
20 
21 
22 //---------------------------------------------------------------------
23 void CRTHitRecoAlg::reconfigure(const Config& config){
24  fVerbose = config.Verbose();
25  fUseReadoutWindow = config.UseReadoutWindow();
26  fQPed = config.QPed();
27  fQSlope = config.QSlope();
28  fPropDelay = config.PropDelay();
29  fPEThresh = config.PEThresh();
30  ftopGain = config.topGain();
31  ftopPed = config.topPed();
32  fSiPMtoFEBdelay = config.SiPMtoFEBdelay();
33  fCoinWindow = config.CoinWindow();
34  fCrtWindow = config.CrtWindow();
35  foutCSVFile = config.outCSVFile();
36  fCSVFile = config.CSVFile();
37  fData = config.Data();
38  if (foutCSVFile)
39  filecsv.open(fCSVFile.c_str());
40  return;
41 }
42 
43 //---------------------------------------------------------------------------------------
44 vector<art::Ptr<CRTData>> CRTHitRecoAlg::PreselectCRTData(const vector<art::Ptr<CRTData>> &crtList, uint64_t trigger_timestamp){
45  if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << "In total " << crtList.size() << " CRTData found in an event" << '\n';
46  vector<art::Ptr<CRTData>> crtdata;
47  bool presel = false;
48 
49  for (size_t febdat_i=0; febdat_i<crtList.size(); febdat_i++) {
50 
51  uint8_t mac = crtList[febdat_i]->fMac5;
52  int adid = fCrtutils->MacToAuxDetID(mac,0);
53  char type = fCrtutils->GetAuxDetType(adid);
54 
55  /// Looking for data within +/- 3ms within trigger time stamp
56  /// Here t0 - trigger time -ve
57  if (fData && (std::fabs(int64_t(crtList[febdat_i]->fTs0 - trigger_timestamp)) > fCrtWindow)) continue;
58 
59  if ( type == 'm'){
60  for(int chan=0; chan<32; chan++) {
61  std::pair<double,double> const chg_cal = fChannelMap->getSideCRTCalibrationMap((int)crtList[febdat_i]->fMac5,chan);
62  float pe = (crtList[febdat_i]->fAdc[chan]-chg_cal.second)/chg_cal.first;
63  if(pe<=fPEThresh) continue;
64  presel = true;
65  if(fVerbose)
66  mf::LogInfo("CRTHitRecoAlg: ") << "\nfebP (mac5, channel, gain, pedestal, adc, pe) = (" << (int)crtList[febdat_i]->fMac5 << ", " << chan << ", "
67  << chg_cal.first << ", " << chg_cal.second << "," << crtList[febdat_i]->fAdc[chan] << "," << pe << ")\n";
68  }
69  }else if ( type == 'c' ) {
70  for(int chan=0; chan<32; chan++) {
71  //float pe = (crtList[febdat_i]->fAdc[chan]-fQPed)/fQSlope;
72  //if(pe<=fPEThresh) continue;
73  presel = true;
74  }
75  }else if ( type == 'd'){
76  for(int chan=0; chan<64; chan++) {
77  float pe = (crtList[febdat_i]->fAdc[chan]-fQPed)/fQSlope;
78  if(pe<=fPEThresh) continue;
79  presel = true;
80  }
81  }
82 
83  if (presel) crtdata.push_back(crtList[febdat_i]);
84  presel = false;
85 
86  }
87  mf::LogInfo("CRTHitRecoAlg:") << "Found " << crtdata.size() << " after preselection "<< '\n';
88  //std::cout<<trigger_timestamp<<std::endl;
89  return crtdata;
90 }
91 
92 //---------------------------------------------------------------------------------------
93 vector<pair<sbn::crt::CRTHit, vector<int>>> CRTHitRecoAlg::CreateCRTHits(vector<art::Ptr<CRTData>> crtList, uint64_t trigger_timestamp) {
94 
95  vector<pair<CRTHit, vector<int>>> returnHits;
96  vector<int> dataIds;
97 
98  uint16_t nMissC = 0, nMissD = 0, nMissM = 0, nHitC = 0, nHitD = 0, nHitM = 0;
99  if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << "Found " << crtList.size() << " FEB events" << '\n';
100 
101  map<string,int> regCounts;
102  std::set<string> regs;
103  map<string,vector<size_t>> sideRegionToIndices;
104 
105  // sort by the time
106  std::sort(crtList.begin(), crtList.end(), compareBytime);
107 
108  //Load Delays map for Top CRT
109  CRT_delay_map const FEB_delay_map = LoadFEBMap();
110  std::vector<std::pair<int,ULong64_t>> CRTReset;
111  ULong64_t TriggerArray[305]={0};
112  for (size_t crtdat_i=0; crtdat_i<crtList.size(); crtdat_i++) {
113  uint8_t mac = crtList[crtdat_i]->fMac5;
114  int adid = fCrtutils->MacToAuxDetID(mac,0);
115  char type = fCrtutils->GetAuxDetType(adid);
116  //For the time being, Only Top CRT delays are loaded, nothing to do for Side CRT yet
117  if (type == 'c' && crtList[crtdat_i]->IsReference_TS1()) {
118  ULong64_t Ts0T1ResetEvent = crtList[crtdat_i]->fTs0 + FEB_delay_map.at((int)mac+73).T0_delay - FEB_delay_map.at((int)mac+73).T1_delay;
119  TriggerArray[(int) mac]=Ts0T1ResetEvent;
120  CRTReset.emplace_back((int) mac,Ts0T1ResetEvent);
121  }
122  }
123  const int trigger_offset= 60; //Average distance between Global Trigger and Trigger_timestamp (ns)
124  ULong64_t GlobalTrigger= trigger_timestamp;
125  if (!CRTReset.empty()) GlobalTrigger = GetMode(CRTReset);
126  //Add average difference between trigger_timestamp and Global trigger
127  else GlobalTrigger=GlobalTrigger-trigger_offset;// In this event, the T1 Reset was probably "vetoed" by the T0 Reset
128  for (int i=0; i<304; i++){
129  if (TriggerArray[i]==0) TriggerArray[i]=GlobalTrigger;
130  }
131  //std::cout<<"Global Trigger "<<GlobalTrigger<<std::endl;
132  //loop over time-ordered CRTData
133  for (size_t febdat_i=0; febdat_i<crtList.size(); febdat_i++) {
134 
135  uint8_t mac = crtList[febdat_i]->fMac5;
136  int adid = fCrtutils->MacToAuxDetID(mac,0); //module ID
137 
138  string region = fCrtutils->GetAuxDetRegion(adid);
139  char type = fCrtutils->GetAuxDetType(adid);
140  CRTHit hit;
141 
142  dataIds.clear();
143 
144  //CERN modules (intramodule coincidence)
145  if ( type == 'c' ) {
146  hit = MakeTopHit(crtList[febdat_i], TriggerArray);
147  if(IsEmptyHit(hit))
148  nMissC++;
149  else {
150  dataIds.push_back(febdat_i);
151  returnHits.push_back(std::make_pair(hit,dataIds));
152  if ((regs.insert(region)).second) regCounts[region] = 1;
153  else regCounts[region]++;
154 
155  nHitC++;
156  }
157  }
158 
159  //DC modules (intramodule coincidence)
160  if ( type == 'd' ) {
161  hit = MakeBottomHit(crtList[febdat_i]);
162  if(IsEmptyHit(hit))
163  nMissD++;
164  else {
165  dataIds.push_back(febdat_i);
166  returnHits.push_back(std::make_pair(hit,dataIds));
167  if ((regs.insert(region)).second) regCounts[region] = 1;
168  else regCounts[region]++;
169 
170  nHitD++;
171  }
172  }
173 
174  if ( type == 'm' )
175  sideRegionToIndices[region].push_back(febdat_i);
176 
177  }//loop over CRTData products
178 
179  vector<size_t> unusedDataIndex;
180  for(auto const& regIndices : sideRegionToIndices) {
181 
182  if(fVerbose)
183  mf::LogInfo("CRTHitRecoAlg: ") << "searching for side CRT hits in region, " << regIndices.first << '\n';
184 
185  vector<size_t> indices = regIndices.second;
186 
187  if(fVerbose)
188  mf::LogInfo("CRTHitRecoAlg: ") << "number of hits associated to this region : " << indices.size() << '\n';
189 
190  for(size_t index_i=0; index_i < indices.size(); index_i++) {
191 
192  dataIds.clear();
193  dataIds.push_back(indices[index_i]);
194  vector<art::Ptr<CRTData>> coinData = {crtList[indices[index_i]]};
195 
196  if(fVerbose)
197  mf::LogInfo("CRTHitRecoAlg: ") << "size .. " << coinData.size()
198  << " data products enetring to time ordring" << '\n';
199 
200  //inner loop over data after data_i in time
201  for (size_t index_j=index_i+1; index_j<indices.size(); index_j++) {
202 
203  if(fVerbose)
204  mf::LogInfo("CRTHitRecoAlg: ") << "i \t"<<index_i << ", j \t"<<index_j << "\t"<<crtList[indices[index_j]]->fTs0 << "\t"<<crtList[indices[index_i]]->fTs0
205  << "\t"<<crtList[indices[index_i]]->fTs0+ fCoinWindow <<'\n';
206 
207  if(crtList[indices[index_j]]->fTs0 < crtList[indices[index_i]]->fTs0)
208  mf::LogError("CRTHitRecoAlg::CreateCRTHits") <<
209  "bad time ordering!" << '\n';
210 
211  if(fVerbose)
212  mf::LogInfo("CRTHitRecoAlg: ") << "size .. " << coinData.size()
213  << " data products before coincidence" << '\n';
214  // in coincidence
215  // if(crtList[indices[index_j]]->fTs0 <= crtList[indices[index_i]]->fTs0 + fCoinWindow) {
216  // if(std::llabs(crtList[indices[index_j]]->fTs0 - crtList[indices[index_i]]->fTs0) < fCoinWindow) {
217  if( (crtList[indices[index_j]]->fTs0>=crtList[indices[index_i]]->fTs0 &&
218  (crtList[indices[index_j]]->fTs0 - crtList[indices[index_i]]->fTs0) < fCoinWindow) ||
219  (crtList[indices[index_j]]->fTs0<crtList[indices[index_i]]->fTs0 &&
220  (crtList[indices[index_i]]->fTs0 - crtList[indices[index_j]]->fTs0) < fCoinWindow)) {
221  if(fVerbose)
222  mf::LogInfo("CRTHitRecoAlg: ") << " in coincidence: i \t " << index_i << " ,j: \t" << index_j <<",i mac: \t"
223  << (int)crtList[indices[index_i]]->fMac5 << ", j mac: \t" <<(int)crtList[indices[index_j]]->fMac5<< '\n';
224 
225  coinData.push_back(crtList[indices[index_j]]);
226  dataIds.push_back(indices[index_j]);
227  }
228 
229  //out of coinWindow
230 
231  if ((crtList[indices[index_j]]->fTs0 - crtList[indices[index_i]]->fTs0) > fCoinWindow || index_j==indices.size()-1){
232  if(fVerbose)
233  mf::LogInfo("CRTHitRecoAlg: ") << "out of coincidence " << index_j << "\t" << indices.size() <<"\t" <<indices.size()-1
234  << " data products..." << '\n';
235  if(fVerbose)
236  mf::LogInfo("CRTHitRecoAlg: ") << "attempting to produce MINOS hit from " << coinData.size()
237  << " data products..." << '\n';
238 
239  CRTHit hit = MakeSideHit(coinData, TriggerArray);
240 
241  if(IsEmptyHit(hit)){
242  unusedDataIndex.push_back(indices[index_i]);
243  nMissM++;
244  }
245  else {
246  if(fVerbose)
247  mf::LogInfo("CRTHitRecoAlg: ") << "MINOS hit produced" << '\n';
248 
249  returnHits.push_back(std::make_pair(hit,dataIds));
250 
251  if ((regs.insert(regIndices.first)).second)
252  regCounts[regIndices.first] = 1;
253  else
254  regCounts[regIndices.first]++;
255 
256  nHitM++;
257  }
258  index_i = index_j-1;
259  if(index_j==indices.size()-1)
260  index_i++;
261 
262  break;
263  }//if jth data out of coinc window
264  }//inner loop over data
265  }// outer loop over data
266  }//loop over side CRTData products
267 
268  if(fVerbose) {
269  mf::LogInfo("CRT") << returnHits.size() << " CRT hits produced!" << '\n'
270  << " nHitC: " << nHitC << " , nHitD: " << nHitD << " , nHitM: " << nHitM << '\n'
271  << " nMisC: " << nMissC << " , nMisD: " << nMissD << " , nMisM: " << nMissM << '\n';
272  auto cts = regCounts.begin();
273  mf::LogInfo("CRT") << " CRT Hits by region" << '\n';
274  while (cts != regCounts.end()) {
275  mf::LogInfo("CRT") << "reg: " << (*cts).first << " , hits: " << (*cts).second << '\n';
276  cts++;
277  }
278  }//if Verbose
279 
280  return returnHits;
281 }
282 //--------------------------------------------------------------------------------------------
283 // Function to make filling a CRTHit a bit faster
284 sbn::crt::CRTHit CRTHitRecoAlg::FillCRTHit(vector<uint8_t> tfeb_id, map<uint8_t,vector<pair<int,float>>> tpesmap,
285  float peshit, uint64_t time0, Long64_t time1, int plane,
286  double x, double ex, double y, double ey, double z, double ez, string tagger){
287  CRTHit crtHit;
288  crtHit.feb_id = tfeb_id;
289  crtHit.pesmap = tpesmap;
290  crtHit.peshit = peshit;
291  crtHit.ts0_s_corr = time0 / 1'000'000'000;
292  crtHit.ts0_ns = time0 % 1'000'000'000;
293  crtHit.ts0_ns_corr = time0;
294  crtHit.ts1_ns = time1 /*% 1'000'000'000*/; //TODO: Update the CRTHit data product /sbnobj/common/CRT . Discussion with SBND people needed
295  crtHit.ts0_s = time0 / 1'000'000'000;
296  crtHit.plane = plane;
297  crtHit.x_pos = x;
298  crtHit.x_err = ex;
299  crtHit.y_pos = y;
300  crtHit.y_err = ey;
301  crtHit.z_pos = z;
302  crtHit.z_err = ez;
303  crtHit.tagger = tagger;
304 
305  return crtHit;
306 
307 } // CRTHitRecoAlg::FillCRTHit()
308 
309 //------------------------------------------------------------------------------------------
310 sbn::crt::CRTHit CRTHitRecoAlg::MakeTopHit(art::Ptr<CRTData> data, ULong64_t GlobalTrigger[305]){
311 
312  uint8_t mac = data->fMac5;
313  if(fCrtutils->MacToType(mac)!='c')
314  mf::LogError("CRTHitRecoAlg::MakeTopHit")
315  << "CRTUtils returned wrong type!" << '\n';
316 
317  map< uint8_t, vector< pair<int,float> > > pesmap;
318  int adid = fCrtutils->MacToAuxDetID(mac,0); //module ID
319  auto const& adGeo = fGeometryService->AuxDet(adid); //module
320  string region = fCrtutils->GetAuxDetRegion(adid);
321  int plane = fCrtutils->AuxDetRegionNameToNum(region);
322  double hitpoint[3], hitpointerr[3], hitlocal[3];
323  TVector3 hitpos (0.,0.,0.);
324  float petot = 0., pemax=0., pemaxx=0., pemaxz=0.;
325  int adsid_max = -1, nabove=0;
326  TVector3 postrig;
327  bool findx = false, findz = false;
328  int maxx=0, maxz=0;
329  double sum=0;
330  for(int chan=0; chan<32; chan++) {
331  sum=sum+data->fAdc[chan];
332  float pe = (data->fAdc[chan]-ftopPed)/ftopGain;
333 // if(pe<=fPEThresh) continue;
334  nabove++;
335  int adsid = fCrtutils->ChannelToAuxDetSensitiveID(mac,chan);
336  petot += pe;
337  pesmap[mac].push_back(std::make_pair(chan,pe));
338 
339  //TVector3 postmp = fCrtutils->ChanToLocalCoords(mac,chan);
340  //strip along z-direction
341  if(adsid >= 0 && adsid < 8){
342  //hitpos.SetX(pe*postmp.X()+hitpos.X());
343  //hitpos.SetX(postmp.X()+hitpos.X());
344  if(pe>pemaxx){
345  pemaxx = pe;
346  maxx = adsid;
347  }
348  findx = true;
349  }
350  //strip along x-direction
351  else if(adsid >= 8 && adsid < 16 ){
352  //hitpos.SetZ(pe*postmp.Z()+hitpos.Z());
353  //hitpos.SetZ(postmp.Z()+hitpos.Z());
354  if(pe > pemaxz) {
355  pemaxz = pe;
356  maxz = adsid;
357  }
358  findz = true;
359  }
360  else {
361  mf::LogError("CRTHitRecoAlg::MakeTopHit")
362  << "auxDetSensitive ID out of range!" << '\n';
363  }
364  //identify trigger channel
365  if(pe>pemax) {
366  TVector3 postmp = fCrtutils->ChanToLocalCoords(mac,chan);
367  adsid_max = chan;
368  pemax = pe;
369  postrig = postmp;
370  }
371  }
372 
373  TVector3 postmp = fCrtutils->ChanToLocalCoords(mac,maxx*2);
374  hitpos.SetX(postmp.X());
375  postmp = fCrtutils->ChanToLocalCoords(mac,maxz*2);
376  hitpos.SetZ(postmp.Z());
377 
378  if(!findx)
379  mf::LogWarning("CRTHitRecoAlg::MakeTopHit") << " no interlayer coincidence found! Missing X coord." << '\n';
380  if(!findz)
381  mf::LogWarning("CRTHitRecoAlg::MakeTopHit") << " no interlayer coincidence found! Missing Z coord." << '\n';
382 
383  //no channels above threshold? return empty hit
384  if(nabove==0||!findx||!findz)
385  return FillCRTHit({},{},0,0,0,0,0,0,0,0,0,0,"");
386 
387  //hitpos*=1.0/petot; //hit position weighted by deposited charge
388  //hitpos*=1.0/nabove;
389  hitlocal[0] = hitpos.X();
390  hitlocal[1] = 0.;
391  hitlocal[2] = hitpos.Z();
392 
393  auto const& adsGeo = adGeo.SensitiveVolume(adsid_max); //trigger strip
394  uint64_t thit = data->fTs0;
395  Long64_t thit1 = data->fTs1;
396 
397  if(adsid_max<8){
398  thit -= (uint64_t)round(abs((92+hitpos.X())*fPropDelay));
399  thit1 -= (uint64_t)round(abs((92+hitpos.X())*fPropDelay));
400  thit -= fSiPMtoFEBdelay; //Correction for 12 ns signal cable from SiPM to FEB
401  thit1 -= fSiPMtoFEBdelay; //Correction for 12 ns signal cable from SiPM to FEB
402  }
403  else{
404  thit -= (uint64_t)round(abs((92+hitpos.Z())*fPropDelay));
405  thit1 -= (uint64_t)round(abs((92+hitpos.Z())*fPropDelay));
406  thit -= fSiPMtoFEBdelay; //Correction for 12 ns signal cable from SiPM to FEB
407  thit1 -= fSiPMtoFEBdelay; //Correction for 12 ns signal cable from SiPM to FEB
408  }
409 
410  adGeo.LocalToWorld(hitlocal,hitpoint); //tranform from module to world coords
411 
412  hitpointerr[0] = adsGeo.HalfWidth1()*2/sqrt(12);
413  hitpointerr[1] = adGeo.HalfHeight();
414  hitpointerr[2] = adsGeo.HalfWidth1()*2/sqrt(12);
415  thit1 = (Long64_t)(thit-GlobalTrigger[(int)mac+73]);
416 
417  //Remove T1 Reset event not correctly flagged, remove T1 reset events, remove T0 reset events
418  if((sum<10000 && thit1<2'001'000 && thit1>2'000'000)||data->IsReference_TS1() || data->IsReference_TS0()) return FillCRTHit({},{},0,0,0,0,0,0,0,0,0,0,"");
419 
420  CRTHit hit = FillCRTHit({mac},pesmap,petot,thit,thit1,plane,hitpoint[0],hitpointerr[0],
421  hitpoint[1],hitpointerr[1],hitpoint[2],hitpointerr[2],region);
422 
423  return hit;
424 
425 } // CRTHitRecoAlg::MakeTopHit
426 
427 //------------------------------------------------------------------------------------------
428 sbn::crt::CRTHit CRTHitRecoAlg::MakeBottomHit(art::Ptr<CRTData> data){
429 
430  uint8_t mac = data->fMac5;
431  map< uint8_t, vector< pair<int,float> > > pesmap;
432  int adid = fCrtutils->MacToAuxDetID(mac,0); //module ID
433  auto const& adGeo = fGeometryService->AuxDet(adid); //module
434  string region = fCrtutils->GetAuxDetRegion(adid);
435  int plane =fCrtutils->AuxDetRegionNameToNum(region);
436  double hitpoint[3], hitpointerr[3], hitlocal[3];
437  TVector3 hitpos (0.,0.,0.);
438  float petot = 0., pemax=0.;
439  int adsid_max = -1, nabove=0;
440  TVector3 postrig;
441  double xmin=0.,xmax=0.;
442 
443  for(int chan=0; chan<64; chan++) {
444 
445  float pe = (data->fAdc[chan]-fQPed)/fQSlope;
446  if(pe<=fPEThresh) continue;
447  nabove++;
448  int adsid = fCrtutils->ChannelToAuxDetSensitiveID(mac,chan);
449  petot += pe;
450  pesmap[mac].push_back(std::make_pair(chan,pe));
451 
452  TVector3 postmp = fCrtutils->ChanToLocalCoords(mac,chan);
453  //all strips along z-direction
454  hitpos.SetX(pe*postmp.X()+hitpos.X());
455  if(postmp.X()<xmin)
456  xmin = postmp.X();
457  if(postmp.X()>xmax)
458  xmax = postmp.X();
459 
460  //identify trigger channel
461  if(pe>pemax) {
462  adsid_max = adsid;
463  pemax = pe;
464  postrig = postmp;
465  }
466  }
467 
468  //no channels above threshold? return empty hit
469  if(nabove==0)
470  return FillCRTHit({},{},0,0,0,0,0,0,0,0,0,0,"");
471 
472  hitpos*=1.0/petot; //hit position weighted by deposited charge
473  hitlocal[0] = hitpos.X();
474  hitlocal[1] = 0.;
475  hitlocal[2] = 0;
476 
477  auto const& adsGeo = adGeo.SensitiveVolume(adsid_max); //trigger strip
478  uint64_t thit = data->fTs0 - adsGeo.HalfLength()*fPropDelay;
479 
480  adGeo.LocalToWorld(hitlocal,hitpoint); //tranform from module to world coords
481 
482  hitpointerr[0] = (xmax-xmin+2*adsGeo.HalfWidth1()*2)/sqrt(12);
483  hitpointerr[1] = adGeo.HalfHeight();
484  hitpointerr[2] = adsGeo.Length()/sqrt(12);
485 
486  CRTHit hit = FillCRTHit({mac},pesmap,petot,thit,thit,plane,hitpoint[0],hitpointerr[0],
487  hitpoint[1],hitpointerr[1],hitpoint[2],hitpointerr[2],region);
488 
489  return hit;
490 
491 } // CRTHitRecoAlg::MakeBottomHit
492 
493 //-----------------------------------------------------------------------------------
494 sbn::crt::CRTHit CRTHitRecoAlg::MakeSideHit(vector<art::Ptr<CRTData>> coinData, ULong64_t GlobalTrigger[305]) {
495 
496  vector<uint8_t> macs;
497  map< uint8_t, vector< pair<int,float> > > pesmap;
498 
499 
500  struct infoA {
501  uint8_t mac5s;
502  int channel;
503  uint64_t t0;
504  TVector3 pos;
505  int strip;
506  };
507 
508  struct infoB {
509  uint8_t mac5s;
510  int channel;
511  uint64_t t0;
512  TVector3 pos;
513  int strip;
514  };
515 
516  vector<infoA> informationA;
517  vector<infoA> informationB;
518 
519  int adid = fCrtutils->MacToAuxDetID(coinData[0]->fMac5,0); //module ID
520  auto const& adGeo = fGeometryService->AuxDet(adid); //module
521  string region = fCrtutils->GetAuxDetRegion(adid);
522  int plane =fCrtutils->AuxDetRegionNameToNum(region);
523 
524 
525  double hitpoint[3], hitpointerr[3];
526  TVector3 hitpos (0.,0.,0.);
527 
528 
529  float petot = 0., pemax = 0., pex=0., pey=0.;
530  int adsid_max = -1, nabove=0, nx=0, ny=0, nz = 0, ntrig = 0;
531 
532  TVector3 postrig;
533 
534  vector<uint64_t> ttrigs;
535  vector<uint64_t> t1trigs;
536  vector<TVector3> tpos;
537  double zmin=DBL_MAX, zmax = -DBL_MAX;
538  double ymin=DBL_MAX, ymax = -DBL_MAX;
539  double xmin=DBL_MAX, xmax = -DBL_MAX;
540  std::vector<int> layID;
541  std::vector<int> febA;
542  std::vector<int> febB;
543 
544 
545  uint64_t southt0_v = -999, southt0_h =-999;
546 
547  //loop over FEBs
548  for(auto const& data : coinData) {
549 
550  if (adid == (int)fCrtutils->MacToAuxDetID((int)data->fMac5,0)){
551  febA.push_back(data->fMac5);
552  }else {
553  febB.push_back(data->fMac5);
554  }
555  }
556 
557  if(fVerbose)
558  std ::cout << "line 451: size of febA: \t" << (int)febA.size()
559  << " size of febB: " << (int)febB.size() << '\n';
560 
561 
562  //loop over FEBs
563  for(auto const& data : coinData) {
564 
565  //if(!(region=="South")) continue;
566  macs.push_back(data->fMac5);
567  adid = fCrtutils->MacToAuxDetID(macs.back(),0);
568 
569 
570  int layer = fCrtutils->GetMINOSLayerID(adid);
571  layID.push_back(layer);
572 
573  auto idx = &data - coinData.data();
574  if(fVerbose) std :: cout << "index: " << idx <<" , mac5: " << (int)macs.back() << " ,adid: " << adid<< '\n';
575 
576 
577  if ((int)febA.size() == 0 or (int)febB.size() == 0) continue;
578 
579  // deciding to use largest size of FEBs in the loop
580  int size = 0;
581  if (febA.size() < febB.size()) size = febB.size();
582  else size = febA.size();
583 
584 
585  for (int i = 0; i < size; i++){
586  if (macs.back() == (int)febA[i]) {
587  //loop over channels
588 
589  for(int chan=0; chan<32; chan++) {
590  std::pair<double,double> const chg_cal = fChannelMap->getSideCRTCalibrationMap((int)data->fMac5,chan);
591  float pe = (data->fAdc[chan]-chg_cal.second)/chg_cal.first;
592 
593  //float pe = (data->fAdc[chan]-fQPed)/fQSlope;
594  if(pe<=fPEThresh) continue;
595  nabove++;
596  if(fVerbose)
597  mf::LogInfo("CRTHitRecoAlg: ") << "\nfebA (mac5, channel, gain, pedestal, adc, pe) = (" << (int)data->fMac5 << ", " << chan << ", "
598  << chg_cal.first << ", " << chg_cal.second << "," << data->fAdc[chan] << "," << pe << ")\n";
599  int adsidA = fCrtutils->ChannelToAuxDetSensitiveID(macs.back(),chan);
600  TVector3 postmp = fCrtutils->ChanToWorldCoords(macs.back(),chan);
601 
602  informationA.push_back ({macs.back(), chan, data->fTs0, postmp, adsidA});
603  }
604 
605  if(fVerbose)
606  mf::LogInfo("CRTHitRecoAlg: ") << "looking for mac " << (int)macs.back()
607  << " and matching with febA : " << (int)febA[i]<<'\n';
608  }else if ( macs.back() == (int)febB[i]){
609 
610  //loop over channels
611  for(int chan=0; chan<32; chan++) {
612  std::pair<double,double> const chg_cal = fChannelMap->getSideCRTCalibrationMap((int)data->fMac5,chan);
613  float pe = (data->fAdc[chan]-chg_cal.second)/chg_cal.first;
614  //float pe = (data->fAdc[chan]-fQPed)/fQSlope;
615  if(pe<=fPEThresh) continue;
616  nabove++;
617  if(fVerbose)
618  mf::LogInfo("CRTHitRecoAlg: ") << "\nfebB (mac5, channel, gain, pedestal, adc, pe) = (" << (int)data->fMac5 << ", " << chan << ", "
619  << chg_cal.first << ", " << chg_cal.second << "," << data->fAdc[chan] << "," << pe << ")\n";
620  int adsidB = fCrtutils->ChannelToAuxDetSensitiveID(macs.back(),chan);
621  TVector3 postmp = fCrtutils->ChanToWorldCoords(macs.back(),chan);
622 
623  informationB.push_back ({macs.back(), chan, data->fTs0, postmp, adsidB});
624  }
625  if(fVerbose)
626  mf::LogInfo("CRTHitRecoAlg: ") << "else if looking for mac "<< (int)macs.back()
627  << " and matching with febB : " << (int)febB[i]<<'\n';
628  }
629  }
630 
631  if(fVerbose)
632  mf::LogInfo("CRTHitRecoAlg: ") << "In CRTHitRecoAlg::MakeSideHit 1st feb is " << (int)fCrtutils->ADToMac(adid).first
633  << " ,2nd feb :"<< (int)fCrtutils->ADToMac(adid).second
634  << ", time "<< data->fTs0
635  << " with module number " << adid << ", no. of FEB " <<'\n';
636 
637  if(fVerbose)
638  mf::LogInfo("CRTHitRecoAlg: ") << "In CRTHitRecoAlg::MakeSideHit functions mac is " << (int)macs.back()
639  << " with module number " << adid << ", no. of FEB " <<'\n';
640 
641  //loop over channels
642  for(int chan=0; chan<32; chan++) {
643  // std :: cout << "chan: ---------------- " << chan << " , "<< fCrtutils->ChannelToAuxDetSensitiveID(macs.back(),chan) << '\n';
644  std::pair<double,double> const chg_cal = fChannelMap->getSideCRTCalibrationMap((int)data->fMac5,chan);
645  float pe = (data->fAdc[chan]-chg_cal.second)/chg_cal.first;
646  //float pe = (data->fAdc[chan]-fQPed)/fQSlope;
647  if(pe<=fPEThresh) continue;
648  nabove++;
649  if(fVerbose)
650  mf::LogInfo("CRTHitRecoAlg: ") << "\nfebC (mac5, channel, gain, pedestal, adc, pe) = (" << (int)data->fMac5 << ", " << chan << ", "
651  << chg_cal.first << ", " << chg_cal.second << "," << data->fAdc[chan] << "," << pe << ")\n";
652  int adsid = fCrtutils->ChannelToAuxDetSensitiveID(macs.back(),chan);
653  petot += pe;
654  pesmap[macs.back()].push_back(std::make_pair(chan,pe));
655 
656  TVector3 postmp = fCrtutils->ChanToWorldCoords(macs.back(),chan);
657 
658  if(fVerbose)
659  mf::LogInfo("CRTHitRecoAlg: ") << "feb: " << (int)macs.back() << " ,chan : \t" << chan
660  << " ,pe: \t"<< pe << ", adc:\t" << data->fAdc[chan] << ", time: \t"<< data->fTs0
661  << " ,x: \t"<< postmp.X() <<" ,y: \t" << postmp.Y() <<" ,z: \t" << postmp.Z() << " petotal : "<< petot << '\n';
662 
663  //East/West Walls (all strips along z-direction) or
664  // North/South inner walls (all strips along x-direction)
665  // All the horizontal layers measure Y first,
666  if(!(region=="South" && layer==1)) {
667  //hitpos.SetY(pe*postmp.Y()+hitpos.Y());
668  //southvertypos.SetX(pe*postmp.X()+southvertypos.X());
669  hitpos.SetY(1.0*postmp.Y()+hitpos.Y());
670  if (fVerbose){
671  mf::LogInfo("CRTHitRecoAlg: ") << "!(region==South && layer==1) : \t" << " feb: " << (int)macs.back() << " ,chan : \t" << chan
672  << " ,pe: \t"<< pe << ", adc:\t" << data->fAdc[chan] << ", time: \t"<< data->fTs0
673  << " ,x: \t"<< postmp.X() <<" ,y: \t" << postmp.Y() <<" ,z: \t" << postmp.Z() << " petotal : "<< petot<< '\n';
674  }
675  ny++;
676  pey+=pe;
677  if(postmp.Y()<ymin)
678  ymin = postmp.Y();
679  if(postmp.Y()>ymax)
680  ymax = postmp.Y();
681  if(region!="South") { //region is E/W/N
682  // hitpos.SetX(pe*postmp.X()+hitpos.X());
683  hitpos.SetX(1.0*postmp.X()+hitpos.X());
684  nx++;
685  pex+=pe;
686  if(postmp.X()<xmin)
687  xmin = postmp.X();
688  if(postmp.X()>xmax)
689  xmax = postmp.X();
690  }
691  }
692  else { //else vertical strips in South wall
693  // hitpos.SetX(pe*postmp.X()+hitpos.X());
694  // hitpos.SetY(pe*postmp.Y()+hitpos.Y());
695  // southvertypos.SetY(pe*postmp.Y()+southvertypos.Y());
696  hitpos.SetX(1.0*postmp.X()+hitpos.X());
697  if (fVerbose){
698  mf::LogInfo("CRTHitRecoAlg: ") << "vertical strips in South wall : \t" << " feb: " << (int)macs.back() << " ,chan : \t" << chan
699  << " ,pe: \t"<< pe << ", adc:\t" << data->fAdc[chan] << ", time: \t"<< data->fTs0
700  << " ,x: \t"<< postmp.X() <<" ,y: \t" << postmp.Y() <<" ,z: \t" << postmp.Z() << " petotal : "<< petot<< '\n';
701  }
702  nx++;
703  pex+=pe;
704  if(postmp.X()<xmin)
705  xmin = postmp.X();
706  if(postmp.X()>xmax)
707  xmax = postmp.X();
708  }
709 
710  //nz = ny
711  //hitpos.SetZ(pe*postmp.Z()+hitpos.Z());
712  hitpos.SetZ(1.0*postmp.Z()+hitpos.Z());
713  nz++;
714  if (fVerbose){
715  if(region =="South")
716  mf::LogInfo("CRTHitRecoAlg: ") << " South wall z: \t" << " feb: " << (int)macs.back() << " ,chan : \t" << chan
717  << " ,pe: \t"<< pe << ", adc:\t" << data->fAdc[chan] << ", time: \t"<< data->fTs0
718  << " ,x: \t"<< postmp.X() <<" ,y: \t" << postmp.Y() <<" ,z: \t" << postmp.Z()<< " petotal : "<< petot<< '\n';
719  }
720  if(postmp.Z()<zmin)
721  zmin = postmp.Z();
722  if(postmp.Z()>zmax)
723  zmax = postmp.Z();
724 
725  //identify trigger channel
726  if(pe>pemax) {
727  adsid_max = adsid;
728  pemax = pe;
729  postrig = postmp;
730  }
731 
732  }//loop over channels
733 
734 
735 
736  //correct trigger time for propegation delay
737  auto const& adsGeo = adGeo.SensitiveVolume(adsid_max); //trigger stripi
738  //auto const& adsGeo = adGeo.SensitiveVolume(last_chan); //trigger strip
739 
740  // Timing determination
741  // ttrigs[layer].push_back(data->fTs0);// - adsGeo.HalfLength()*fPropDelay);
742 
743 
744  //ttrigs.push_back(data->fTs0 - uint64_t(adsGeo.HalfLength()*fPropDelay));
745  //t1trigs.push_back(data->fTs1 - uint64_t(adsGeo.HalfLength()*fPropDelay));
746  ttrigs.push_back(data->fTs0);
747  t1trigs.push_back(data->fTs1);
748  tpos.push_back(postrig);
749  ntrig++;
750 
751  if(fVerbose)
752  mf::LogInfo("CRTHitRecoAlg: ") << "raw T0: " << data->fTs0 << " ,half lenth : \t" << adsGeo.HalfLength() << " ,delay: \t" <<fPropDelay
753  << " ,corrected time: " << data->fTs0 - uint64_t(adsGeo.HalfLength()*fPropDelay) << '\n';
754 
755 
756  if(region=="South" && layer==1) {
757  southt0_h = data->fTs0;
758  if(fVerbose)
759  mf::LogInfo("CRTHitRecoAlg: ") << "southt0_h : " << layer << "\t" << southt0_h << '\n';
760  }else if (region=="South" && layer!=1){
761  southt0_v = data->fTs0;
762  if(fVerbose)
763  mf::LogInfo("CRTHitRecoAlg: ") << "southt0_v : " << layer << "\t"<< southt0_v << '\n';
764  }else {
765  southt0_h = -999;
766  southt0_v = -999;
767  }
768 
769 
770  }//loop over FEBs
771 
772 
773 
774  std::sort(layID.begin(), layID.end());
775  layID.resize(std::distance(layID.begin(), std::unique(layID.begin(), layID.end())));
776 
777  if(fVerbose)
778  mf::LogInfo("CRTHitRecoAlg: ") << "size of layer ID : \t" <<layID.size() << '\n';
779 
780 
781 
782  //no channels above threshold or no intermodule coincidences? return empty hit
783  // if(nabove==0 || layID.size()!=2) {
784  if(nabove==0 || layID.size() < 2) {
785  if(nabove==0 && fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << "no channels above threshold!" << '\n';
786  if(layID.size()<2 && fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << "no coincidence found" << '\n';
787  return FillCRTHit({},{},0,0,0,0,0,0,0,0,0,0,"");
788  }
789 
790 
791  //-----------------------------------------------------------
792  // Matching between both end readouts for full length module
793  //-----------------------------------------------------------
794  //1. The axis you are measuring on, with the correct sign (d)
795  //2. The centre of the detector (c)
796  //3. The absolute position is c + d · (x - L/2),
797  // that is a displacement from the centre c by a distance x - L/2 (relative to the centre position L/2).
798  // we are determining along the z which is here d = z (aligned with the long axis of the ICARUS detector).
799  // The centre of the module you are reading is adsGeo.GetCenter()
800  // so adsGeo.GetCenter() + geo::Zaxis() * (x - adsGeo.HalfLength()) may work.
801  // We can define a different x, from the middle of the detector:
802 
803  // T1 = t_hit + (L/2 + x)/v_propagation;
804  // T2 = t_hit + (L/2 - x)/v_propagation;
805  // double deltat(T1-T2) = 2*x/v_propagation;
806  //---------------------------------------------------------------------
807  //std::ofstream o(filename.c_str());
808 
809  geo::Point_t posA; geo::Point_t posB;
810  bool layer1 = false;
811  bool layer2 = false;
812 
813  geo::Point_t crossfebpos;
814 
815  uint64_t t0_1=-5; uint64_t t0_2=-5;
816  uint64_t t1_1=-5; uint64_t t1_2=-5;
817  uint64_t t2_1=-5; uint64_t t2_2=-5;
818  int mac5_1=-5; int mac5_2 = -5;
819  geo::Point_t center;
820 
821  for(auto const& infn: informationA) {
822  auto i = &infn - informationA.data();
823  auto const& adsGeo = adGeo.SensitiveVolume(infn.strip); //trigger stripi
824  center = adsGeo.GetCenter();
825 
826  if(fVerbose)
827  mf::LogInfo("CRTHitRecoAlg: ")<< "A type ----------> time: " << (long long int)infn.t0 << " ,macs "<< (int)infn.mac5s //<< '\n';
828  << " ,chal "<< infn.channel
829  << " , position "<<infn.pos[2] << '\n';
830 
831  if ((int)infn.mac5s != (int)informationA[i+1].mac5s and i < (int)informationA.size()-1){
832  layer1 = true;
833 
834  if ((int)infn.mac5s % 2 == 0) t1_1 = infn.t0;
835  else t1_1 = informationA[i+1].t0;
836  if ((int)informationA[i+1].mac5s % 2 != 0) t1_2 = informationA[i+1].t0;
837  else t1_2 = infn.t0;
838  if(fVerbose)
839  mf::LogInfo("CRTHitRecoAlg: ")<< "t1: " << t1_1 << ", t2:"<< t1_2 << ", deltat : "<< int64_t(t1_1 - t1_2) << '\n';
840 
841  float zaxixpos = 0.5*(int64_t(t1_1 - t1_2)/fPropDelay);
842 
843  posA = adsGeo.GetCenter() + geo::Zaxis() * zaxixpos;
844  if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ")<< "posA (==0): "<< posA<< '\n';
845 
846 
847  if(fVerbose)
848  mf::LogInfo("CRTHitRecoAlg: ")<< i << " ,1st mac5: "<< (int)infn.mac5s << " 1st time: " << (long long int)infn.t0
849  << " ,2nd mac5: "<<(int)informationA[i+1].mac5s << ", 2nd time " << (long long int)informationA[i+1].t0 << " , deltaT: "
850  << int64_t(t1_1 - t1_2) << " , length: " << adsGeo.Length()
851  << " ,propagation delay: "<< fPropDelay << " , pos z: " << zaxixpos
852  << " , center: " << adsGeo.GetCenter() << " , zaxis: "<< geo::Zaxis() << " , half length: " << adsGeo.HalfLength()
853  << " , actual pos w.rt. z: " << posA << '\n';
854  }
855  }
856 
857  for(auto const& infn: informationB) {
858 
859 
860  if(fVerbose)
861  mf::LogInfo("CRTHitRecoAlg: ")<< " B type ----------> time: " << infn.t0 << " ,macs "<< (int)infn.mac5s //<< '\n';
862  << " ,chal "<< infn.channel
863  << " , position "<<infn.pos[2] << '\n';
864 
865  auto i = &infn - informationB.data();
866  auto const& adsGeo = adGeo.SensitiveVolume(infn.strip); //trigger stripi
867 
868 
869  if ((int)infn.mac5s != (int)informationB[i+1].mac5s and i < (int)informationB.size()-1){
870 
871  layer2 = true;
872  mac5_2 = (int)infn.mac5s;
873  t0_2 = (uint64_t)infn.t0;
874 
875  if ((int)infn.mac5s % 2 == 0) t2_1 = infn.t0;
876  else t2_1 = informationB[i+1].t0;
877  if ((int)informationB[i+1].mac5s % 2 != 0) t2_2 = informationB[i+1].t0;
878  else t2_2 = infn.t0;
879 
880  if(fVerbose)
881  mf::LogInfo("CRTHitRecoAlg: ")<< "t1: " << t2_1 <<", t2:"<< t2_2 << ", deltat : "<< int64_t(t2_1 - t2_2) << '\n';
882  //if (foutCSVFile) filecsv << plane << "\t"<< int64_t(t2_1 - t2_2) << "\n";
883  float zaxixpos = 0.5*(int64_t(t2_1 - t2_2)/fPropDelay);
884 
885  posB = adsGeo.GetCenter() + geo::Zaxis() * zaxixpos;
886  if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ")<< "posB (== 0): "<< posB<< '\n';
887 
888 
889  if(fVerbose)
890  mf::LogInfo("CRTHitRecoAlg: ")<< i << " ,1st mac5: "<< (int)infn.mac5s << " 1st time: " << (long long int)infn.t0
891  << " ,2nd mac5: "<<(int)informationB[i+1].mac5s << ", 2nd time " << (long long int)informationB[i+1].t0 << " , deltaT: "
892  << int64_t(t2_1 - t2_2) << " , length: " << adsGeo.Length()
893  << " ,propagation delay: "<< fPropDelay << " , pos z: " << zaxixpos
894  << " , center: " << adsGeo.GetCenter() << " , zaxis: "<< geo::Zaxis() << " , half length: " << adsGeo.HalfLength()
895  << " , actual pos w.rt. z: " << posB << '\n';
896 
897  }
898 
899  }
900 
901  int crossfeb = std::abs(mac5_1 - mac5_2);
902 
903 
904  // side crt and match the both layers
905  if (layer1 && layer2 && region!="South" && region!="North" ){//&& nx==4){
906  float avg = 0.5*(posA.Z() + posB.Z());
907  hitpos.SetZ(avg);
908  hitpos.SetX(hitpos.X()*1.0/nx);
909  hitpos.SetY(hitpos.Y()*1.0/nx);
910  //hitpos.SetX(hitpos.X()*1.0/petot);
911  //hitpos.SetY(hitpos.Y()*1.0/petot);
912  if(fVerbose)
913  mf::LogInfo("CRTHitRecoAlg: ") << "z position in layer 1: "<< posA.Z() << " and in layer 2 "<< posB.Z()
914  << " average is "<< (posA.Z()+ posB.Z())/2. << " ,hitpos z " << hitpos[2] << '\n';
915 
916 
917  }else if ((int)informationA.size()==1 and (int)informationB.size()==1
918  and (crossfeb == 7 or crossfeb == 5) and
919  region!="South" && region!="North"){
920 
921  int z_pos = int64_t(t0_1 - t0_2)/(uint64_t(2*fPropDelay));
922  crossfebpos = center + geo::Zaxis()*z_pos;
923 
924  hitpos.SetZ(crossfebpos.Z());
925  hitpos.SetX(hitpos.X()*1.0/nx);
926  hitpos.SetY(hitpos.Y()*1.0/nx);
927  //hitpos.SetX(hitpos.X()*1.0/petot);
928  //hitpos.SetY(hitpos.Y()*1.0/petot);
929  if(fVerbose)
930  mf::LogInfo("CRTHitRecoAlg: ") << "hello hi namaskar, hitpos z " << hitpos[2] << '\n';
931  // side crt and only single layer match
932  }else if (layer1 && region!="South" && region!="North"){// && nx==1){
933  hitpos.SetZ(posA.Z());
934  hitpos.SetX(hitpos.X()*1.0/nx);
935  hitpos.SetY(hitpos.Y()*1.0/nx);
936  //hitpos.SetX(hitpos.X()*1.0/petot);
937  //hitpos.SetY(hitpos.Y()*1.0/petot);
938  if(fVerbose)
939  mf::LogInfo("CRTHitRecoAlg: ") << " same layer coincidence: z position in layer 1: "<< posA.Z() << " ,hitpos z " << hitpos[2] << '\n';
940 
941  // side crt and only single layer match
942  }else if (layer2 && region!="South" && region!="North" ){//&& nx==1){
943  hitpos.SetZ(posB.Z());
944  hitpos.SetX(hitpos.X()*1.0/nx);
945  hitpos.SetY(hitpos.Y()*1.0/nx);
946  //hitpos.SetX(hitpos.X()*1.0/petot);
947  //hitpos.SetY(hitpos.Y()*1.0/petot);
948  if(fVerbose)
949  mf::LogInfo("CRTHitRecoAlg: ") << " same layer coincidence: z position in layer 2 "<< posB.Z() << " ,hitpos z " << hitpos[2] << '\n';
950 
951  }else if (region!="South" && region!="North" ){//&& nx==2){
952  hitpos*=1.0/nx;
953  //hitpos.SetX(hitpos.X()*1.0/petot);
954  //hitpos.SetY(hitpos.Y()*1.0/petot);
955  //hitpos.SetZ(hitpos.Z()*1.0/petot);
956  if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << " In side CRTs [E/W] x: \t"<< hitpos[0] <<" ,y: \t" << hitpos[1] <<" ,z: \t" << hitpos[2]<< '\n';
957  }/*else {
958  hitpos.SetX(-99999);
959  hitpos.SetY(-99999);
960  hitpos.SetZ(-99999);
961  }*/
962 
963 
964  //finish averaging and fill hit point array
965  if(region=="South") {
966  /*
967  hitpos.SetX(hitpos.X()*1.0/pex);
968  hitpos.SetZ(hitpos.Z()*1.0/petot);
969  */
970  ///*
971  hitpos.SetX(hitpos.X()/nx);
972  hitpos.SetY(hitpos.Y()/ny);
973  hitpos.SetZ(hitpos.Z()/nz);
974  // */
975  // }else
976  //hitpos*=1.0/petot; //hit position weighted by deposited charge
977 
978  }else if (region=="North"){
979  //hitpos*=1.0/petot;
980  hitpos*=1.0/nz;
981 
982  //}else if (region!="South" && region!="North"){
983  // hitpos.SetX(hitpos.X()*1.0/petot);
984  //hitpos.SetY(hitpos.Y()*1.0/petot);
985  }
986 
987  //else
988  //hitpos*=1.0/petot; //hit position weighted by deposited charge
989  //hitpos*=1.0; //hit position weighted by deposited charge
990 
991  hitpoint[0] = hitpos.X();
992  hitpoint[1] = hitpos.Y();
993  hitpoint[2] = hitpos.Z();
994 
995  if (region=="South" && hitpoint[0] >= 366. && hitpoint[1] > 200. && fVerbose)
996  mf::LogInfo("CRTHitRecoAlg: ") << "I am looking for south wall : macs " << (int)macs.back() << " x: \t"<< hitpoint[0]
997  <<" ,y: \t" << hitpoint[1] <<" ,z: \t" << hitpoint[2] << '\n';
998 
999  if (fVerbose){
1000  if (region=="North") mf::LogInfo("CRTHitRecoAlg: ") << "north wall x: \t"<< hitpoint[0] <<" ,y: \t" << hitpoint[1] <<" ,z: \t" << hitpoint[2]<< '\n';
1001  }
1002  if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << " nx: \t"<< nx <<" ,ny: \t" << ny <<" ,nz: \t" << nz<< '\n';
1003  if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << " x: \t"<< hitpoint[0] <<" ,y: \t" << hitpoint[1] <<" ,z: \t" << hitpoint[2]<< '\n';
1004 
1005  //time stamp averaged over all FEBs
1006  uint64_t thit = 0., t1hit = 0.;// thit_1 = 0.;
1007 
1008  ttrigs.resize(std::distance(ttrigs.begin(), std::unique(ttrigs.begin(), ttrigs.end())));
1009  //thit = std::reduce(ttrigs.begin(), ttrigs.end()) / (uint64_t)ttrigs.size();
1010  //thit = std::accumulate(ttrigs.begin(), ttrigs.end()) / (uint64_t)ttrigs.size();
1011 
1012  t1trigs.resize(std::distance(t1trigs.begin(), std::unique(t1trigs.begin(), t1trigs.end())));
1013  // t1hit = std::reduce(t1trigs.begin(), t1trigs.end()) / (uint64_t)t1trigs.size();
1014  // t1hit = std::accumulate(t1trigs.begin(), t1trigs.end()) / (uint64_t)t1trigs.size();
1015 
1016  for(uint64_t const t : ttrigs){
1017  if (region=="North" || region=="South") /*thit += t;*/ thit += t-uint64_t(200.*fPropDelay)-fSiPMtoFEBdelay;
1018  else /*thit += t;*/thit += t-uint64_t(400.*fPropDelay)-fSiPMtoFEBdelay;
1019  //thit += t;
1020  //thit += t-uint64_t(400.*fPropDelay);
1021  if (fVerbose)
1022  mf::LogInfo("CRTHitRecoAlg: ") << " Inside the loop: t: \t"<< t << "\t" << uint64_t(200.*fPropDelay)
1023  << " t0hit : "<< thit <<" size ttrig: \t" << ttrigs.size()<< '\n';
1024  }
1025 
1026  thit=thit/uint64_t(ttrigs.size());
1027 
1028  for(double const t1 : t1trigs)
1029  t1hit += t1-uint64_t(400.*fPropDelay)-fSiPMtoFEBdelay;
1030 
1031 
1032  t1hit=t1hit/uint64_t(t1trigs.size());
1033 
1034  if (region=="South" && fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << "..................... Hello ....Welcome to Beam............"<< '\n';
1035  if (fVerbose) mf::LogInfo("CRTHitRecoAlg: ") << " <time>: T0: \t"<< thit << " T1 : "<< t1hit <<" size ttrig: \t" << ttrigs.size()<< '\n';
1036 
1037  if(region=="South" && fVerbose)
1038  mf::LogInfo("CRTHitRecoAlg: ") << "southt0_h: " << southt0_h << " ,southt0_v : " << southt0_v
1039  << " ,deltaT: \t" << int64_t(southt0_h - southt0_v) << '\n';
1040  if (foutCSVFile) filecsv << southt0_h << "\t" << southt0_v << "\t"<< int64_t(southt0_h - southt0_v) << "\n";
1041 
1042  //error estimates (likely need to be revisted)
1043  auto const& adsGeo = adGeo.SensitiveVolume(adsid_max);
1044  if(region!="North" && region!="South"){
1045  hitpointerr[0] = (xmax-xmin)/sqrt(12);
1046  hitpointerr[1] = (ymax-ymin)/sqrt(12);
1047  hitpointerr[2] = (zmax-zmin)/sqrt(12);
1048  // hitpointerr[2] = adsGeo.Length()/sqrt(12);
1049  }
1050 
1051  if(region=="North"){
1052  hitpointerr[0] = (xmax-xmin)/sqrt(12);
1053  hitpointerr[1] = (ymax-ymin)/sqrt(12);
1054  hitpointerr[2] = (zmax-zmin)/sqrt(12);
1055  }
1056 
1057  if(region=="South"){
1058  hitpointerr[0] = adsGeo.HalfWidth1()*2/sqrt(12);
1059  hitpointerr[1] = adsGeo.HalfWidth1()*2/sqrt(12);
1060  hitpointerr[2] = (zmax-zmin)/sqrt(12);
1061  }
1062 
1063  Long64_t thit1=(Long64_t)(thit-GlobalTrigger[(int)macs.at(0)]);
1064 
1065  //generate hit
1066  CRTHit hit = FillCRTHit(macs,pesmap,petot,thit,thit1,plane,hitpoint[0],hitpointerr[0],
1067  hitpoint[1],hitpointerr[1],hitpoint[2],hitpointerr[2],region);
1068 
1069  return hit;
1070 }
1071 
1072 //-----------------------------------------------------------------------------
1073 bool CRTHitRecoAlg::IsEmptyHit(CRTHit hit) {
1074 
1075  if ( hit.feb_id.empty() && hit.pesmap.empty() && hit.peshit == 0
1076  && hit.ts0_ns == 0 && hit.ts1_ns == 0 && hit.plane == 0
1077  && hit.x_pos == 0 && hit.x_err == 0 && hit.y_pos == 0
1078  && hit.y_err == 0 && hit.z_pos == 0 && hit.z_err == 0 && hit.tagger == "")
1079  return true;
1080 
1081  return false;
1082 }
1083 
1084 //-----------------------------------------------------------------------------
1085 ULong64_t icarus::crt::GetMode(std::vector<std::pair<int, ULong64_t>> vector) {
1086 
1087  sort(vector.begin(), vector.end(), icarus::crt::sortbytime);
1088 
1089  int modecounter = 0;
1090  int isnewmodecounter = 0;
1091  ULong64_t Mode = 0;
1092  ULong64_t isnewMode = 0;
1093  bool isFirst = true;
1094  for (auto i : vector) {
1095  if (!isFirst) {
1096  if (i.second == Mode) modecounter++;
1097  else if (i.second !=isnewMode) {
1098  isnewMode = i.second;
1099  isnewmodecounter = 1;
1100  }
1101  else if (i.second == isnewMode) {
1102  isnewmodecounter++;
1103  if (isnewmodecounter > modecounter) {
1104  Mode = isnewMode;
1105  modecounter = isnewmodecounter;
1106  }
1107  }
1108  }
1109  else {
1110  isFirst = false;
1111  Mode = i.second;
1112  modecounter++;
1113  }
1114  }
1115  return Mode;
1116 }
z
process_name opflash particleana ie ie ie z
Definition: run_opflash_electron.fcl:92
icarus::crt::CRTHitRecoAlg::fQSlope
double fQSlope
Pedestal slope of SiPMs [ADC/photon].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:176
sbn::crt::CRTHit::pesmap
std::map< uint8_t, std::vector< std::pair< int, float > > > pesmap
Saves signal hit information (FEB, local-channel and PE) .
Definition: CRTHit.hh:26
icarus::crt::CRTHitRecoAlg::FillCRTHit
CRTHit FillCRTHit(vector< uint8_t > tfeb_id, map< uint8_t, vector< pair< int, float >>> tpesmap, float peshit, uint64_t time0, Long64_t time1, int plane, double x, double ex, double y, double ey, double z, double ez, string tagger)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.cc:284
second
second
Definition: detectorproperties_icarus.fcl:8
ServiceUtil.h
Utilities related to art service access.
sbn::crt::CRTHit::feb_id
std::vector< uint8_t > feb_id
FEB address.
Definition: CRTHit.hh:25
sbn::crt::CRTHit::ts1_ns
double ts1_ns
Timestamp T1 ([signal time w.r.t. Trigger time]), in UTC absolute time scale in nanoseconds from the ...
Definition: CRTHit.hh:34
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
sbn::crt::CRTHit::peshit
float peshit
Total photo-electron (PE) in a crt hit.
Definition: CRTHit.hh:27
icarus::crt::CRTHitRecoAlg::fCSVFile
std::string fCSVFile
FCL input: file name for output CSV File.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:173
icarus::crt::CRTHitRecoAlg::fPEThresh
double fPEThresh
threshold[PE] above which charge amplitudes used in hit reco
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:178
hit
process_name hit
Definition: cheaterreco.fcl:51
icarus::crt::CRTHitRecoAlg::compareBytime
static bool compareBytime(art::Ptr< CRTData > const &a, art::Ptr< CRTData > const &b)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:197
sbn::crt::CRTHit::ts0_ns_corr
double ts0_ns_corr
[Honestly, not sure at this point, it was there since long time (BB)]
Definition: CRTHit.hh:33
icarus::crt::CRTHitRecoAlg::fSiPMtoFEBdelay
uint64_t fSiPMtoFEBdelay
SiPM to FEB cable induced delay: 11.6 [ns].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:181
sbn::crt::CRTHit::ts0_s_corr
double ts0_s_corr
[Honestly, not sure at this point, it was there since long time (BB)]
Definition: CRTHit.hh:30
sbn::crt::CRTHit::ts0_s
uint64_t ts0_s
Second-only part of timestamp T0.
Definition: CRTHit.hh:29
icarus::crt::CRTHitRecoAlg::PreselectCRTData
vector< art::Ptr< CRTData > > PreselectCRTData(const vector< art::Ptr< CRTData >> &crtList, uint64_t trigger_timestamp)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.cc:44
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
icarus::crt::CRT_delay_map
std::map< feb_index, FEB_delay > CRT_delay_map
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:215
icarus::crt::CRTHitRecoAlg::fCoinWindow
uint64_t fCoinWindow
Coincidence window used for grouping side CRT triggers [ns].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:182
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
icarus::crt::CRTHitRecoAlg::MakeTopHit
CRTHit MakeTopHit(art::Ptr< CRTData > data, ULong64_t GlobalTrigger[])
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.cc:310
icarusDB::IICARUSChannelMap::getSideCRTCalibrationMap
virtual std::pair< double, double > getSideCRTCalibrationMap(int mac5, int chan) const =0
geo::vect::indices
constexpr std::array< std::size_t, geo::vect::dimension< Vector >)> indices()
Returns a sequence of indices valid for a vector of the specified type.
Definition: geo_vectors_utils.h:2022
icarus::crt::CRTHitRecoAlg::fUseReadoutWindow
bool fUseReadoutWindow
Only reconstruct hits within TPC readout window.
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:174
icarus::crt::CRTHitRecoAlg::fCrtWindow
uint64_t fCrtWindow
Looking data window within trigger timestamp [ns].
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.h:183
detinfo::trigger_offset
int trigger_offset(DetectorClocksData const &data)
Definition: DetectorClocksData.h:632
icarus::crt::CRTHitRecoAlg::MakeSideHit
CRTHit MakeSideHit(vector< art::Ptr< CRTData >> coinData, ULong64_t GlobalTrigger[])
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.cc:494
crt
process_name crt
Definition: run_genericCRT.fcl:82
icarus::crt::CRTHitRecoAlg::CreateCRTHits
vector< pair< CRTHit, vector< int > > > CreateCRTHits(vector< art::Ptr< CRTData >> crtList, uint64_t trigger_timestamp)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.cc:93
icarus::crt::GetMode
ULong64_t GetMode(std::vector< std::pair< int, ULong64_t >> vector)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTHitRecoAlg.cc:1085
Generated by   doxygen 1.8.5