10 #include "art/Framework/Core/EDAnalyzer.h"
11 #include "art/Framework/Core/ModuleMacros.h"
12 #include "art/Framework/Principal/Event.h"
13 #include "art/Framework/Principal/Handle.h"
14 #include "art/Framework/Principal/Run.h"
15 #include "art/Framework/Principal/SubRun.h"
16 #include "canvas/Utilities/InputTag.h"
17 #include "fhiclcpp/ParameterSet.h"
18 #include "messagefacility/MessageLogger/MessageLogger.h"
19 #include "art_root_io/TFileService.h"
20 #include "art/Framework/Services/Registry/ServiceHandle.h"
21 #include "canvas/Persistency/Common/FindManyP.h"
31 #include "nusimdata/SimulationBase/MCGeneratorInfo.h"
32 #include "nusimdata/SimulationBase/MCTruth.h"
33 #include "nusimdata/SimulationBase/MCParticle.h"
78 void analyze(art::Event
const&
e)
override;
84 bool HitCompare(
const art::Ptr<CRTHit>& h1,
const art::Ptr<CRTHit>& h2);
156 ,fOpHitModuleLabel(
p.get<art::InputTag>(
"OpHitModuleLabel",
"ophit"))
157 ,fOpFlashModuleLabel0(
p.get<art::InputTag>(
"OpFlashModuleLabel0",
"opflashTPC0"))
158 ,fOpFlashModuleLabel1(
p.get<art::InputTag>(
"OpFlashModuleLabel1",
"opflashTPC1"))
161 ,fCrtHitModuleLabel(
p.get<art::InputTag>(
"CrtHitModuleLabel",
"crthit"))
162 ,fTriggerLabel(
p.get<art::InputTag>(
"TriggerLabel",
"daqTrigger"))
164 ,fCoinWindow(
p.get<
double>(
"CoincidenceWindow",60.0))
165 ,fOpDelay(
p.get<
double>(
"OpDelay",55.1))
166 ,fCrtDelay(
p.get<
double>(
"CrtDelay",1.6e6))
167 ,fFlashPeThresh(
p.get<
int>(
"FlashPeThresh",9000))
168 ,fHitPeThresh(
p.get<
int>(
"HitPeThresh",700))
169 ,fFlashVelocity(
p.get<
double>(
"FlashVelocityThresh",-40.))
170 ,fFlashZOffset(
p.get<
double>(
"FlashZOffset",0.))
171 ,fHitVelocityMax(
p.get<
double>(
"HitVelocityMax",20.))
172 ,fHitVelocityMin(
p.get<
double>(
"HitVelocityMin",1.))
178 fFlashLabels[0] = fOpFlashModuleLabel0;
179 fFlashLabels[1] = fOpFlashModuleLabel1;
180 fFlashLabels[2] = fOpFlashModuleLabel2;
181 fFlashLabels[3] = fOpFlashModuleLabel3;
184 fGeometryService = lar::providerFrom<geo::Geometry>();
186 art::ServiceHandle<art::TFileService>
tfs;
188 fMatchTree = tfs->make<TTree>(
"matchTree",
"CRTHit - OpHit/Flash matching analysis");
190 fMatchTree->Branch(
"event", &fEvent,
"event/I");
191 fMatchTree->Branch(
"run", &fRun,
"run/I");
192 fMatchTree->Branch(
"subrun", &fSubRun,
"subrun/I");
193 fMatchTree->Branch(
"ncrt", &fNCrt,
"nCrtHit/I");
194 fMatchTree->Branch(
"crtXYZT", &fCrtXYZT);
195 fMatchTree->Branch(
"crtXYZErr", &fCrtXYZErr);
196 fMatchTree->Branch(
"crtPE", &fCrtPE);
197 fMatchTree->Branch(
"crtRegion", &fCrtRegion);
198 fMatchTree->Branch(
"tofHit", &fTofHit);
199 fMatchTree->Branch(
"tofFlash", &fTofFlash);
200 fMatchTree->Branch(
"tofFlashHit", &fTofFlashHit);
201 fMatchTree->Branch(
"peHit", &fTofPeHit);
202 fMatchTree->Branch(
"peFlash", &fTofPeFlash);
203 fMatchTree->Branch(
"peFlashHit", &fTofPeFlashHit);
204 fMatchTree->Branch(
"xyztHit", &fTofXYZTHit);
205 fMatchTree->Branch(
"xyztFlash", &fTofXYZTFlash);
206 fMatchTree->Branch(
"xyztFlashHit", &fTofXYZTFlashHit);
207 fMatchTree->Branch(
"distHit", &fDistHit);
208 fMatchTree->Branch(
"distFlash", &fDistFlash);
209 fMatchTree->Branch(
"distFlashHit", &fDistFlashHit);
210 fMatchTree->Branch(
"tpcHit", &fTofTpcHit);
211 fMatchTree->Branch(
"tpcFlash", &fTofTpcFlash);
212 fMatchTree->Branch(
"gate_type", &m_gate_type,
"gate_type/b");
213 fMatchTree->Branch(
"gate_name", &m_gate_name);
214 fMatchTree->Branch(
"trigger_timestamp", &m_trigger_timestamp,
"trigger_timestamp/l");
215 fMatchTree->Branch(
"gate_start_timestamp", &m_gate_start_timestamp,
"gate_start_timestamp/l");
216 fMatchTree->Branch(
"trigger_gate_diff", &m_trigger_gate_diff,
"trigger_gate_diff/l");
217 fMatchTree->Branch(
"gate_crt_diff",&m_gate_crt_diff,
"gate_crt_diff/l");
231 MF_LOG_DEBUG(
"CRTPMTMatching: ") <<
"beginning analyis" <<
'\n';
234 fEvent = e.id().event();
236 fSubRun = e.subRun();
241 if( !fTriggerLabel.empty() ) {
243 art::Handle<sbn::ExtraTriggerInfo> trigger_handle;
244 e.getByLabel( fTriggerLabel, trigger_handle );
245 if( trigger_handle.isValid() ) {
247 m_gate_type = (
unsigned int)bit;
249 m_trigger_timestamp = trigger_handle->triggerTimestamp;
250 m_gate_start_timestamp = trigger_handle->beamGateTimestamp;
251 m_trigger_gate_diff = trigger_handle->triggerTimestamp - trigger_handle->beamGateTimestamp;
255 mf::LogError(
"CRTPMTMatching:") <<
"No raw::Trigger associated to label: " << fTriggerLabel.label() <<
"\n" ;
259 std::cout <<
"Trigger Data product " << fTriggerLabel.label() <<
" not found!\n" ;
264 art::Handle< std::vector<recob::OpHit> > opHitListHandle;
265 std::vector< art::Ptr<recob::OpHit> > opHitList;
266 if( e.getByLabel(fOpHitModuleLabel,opHitListHandle) )
267 art::fill_ptr_vector(opHitList, opHitListHandle);
271 map<int, art::Handle< std::vector<recob::OpFlash> > > flashHandles;
272 std::map<int,std::vector< art::Ptr<recob::OpFlash> >> opFlashLists;
274 for(
int i=0; i<2; i++) {
275 if( e.getByLabel(fFlashLabels[i],flashHandles[i]) )
276 art::fill_ptr_vector(opFlashLists[i], flashHandles[i]);
281 art::Handle< std::vector<CRTHit> > crtHitListHandle;
282 std::vector< art::Ptr<CRTHit> > crtHitList;
283 if( e.getByLabel(fCrtHitModuleLabel,crtHitListHandle))
284 art::fill_ptr_vector(crtHitList, crtHitListHandle);
286 fNCrt = crtHitList.size();
288 for(
auto const&
crt : crtHitList){
289 vector<double> xyzt, xyzerr;
290 TVector3 rcrt(
crt->x_pos,
crt->y_pos,
crt->z_pos);
291 m_gate_crt_diff = m_gate_start_timestamp -
crt->ts0_ns;
292 xyzt.push_back(rcrt.X());
293 xyzt.push_back(rcrt.Y());
294 xyzt.push_back(rcrt.Z());
296 double tcrt = double(m_gate_start_timestamp -
crt->ts0_ns)/1e3;
303 xyzt.push_back(tcrt);
304 fCrtXYZT.push_back(xyzt);
306 xyzerr.push_back(
crt->x_err);
307 xyzerr.push_back(
crt->y_err);
308 xyzerr.push_back(
crt->z_err);
309 fCrtXYZErr.push_back(xyzerr);
311 fCrtPE.push_back(
crt->peshit);
312 fCrtRegion.push_back(crtutil->AuxDetRegionNameToNum(
crt->tagger));
316 double tdiff = DBL_MAX, rdiff=DBL_MAX, peflash=DBL_MAX;
319 double flashHitT = DBL_MAX, flashHitPE=DBL_MAX, flashHitDiff=DBL_MAX;
320 vector<double> flashHitxyzt;
322 for(
auto const& flashList : opFlashLists) {
324 art::FindManyP<recob::OpHit> findManyHits(flashHandles[flashList.first], e, fFlashLabels[flashList.first]);
326 for(
size_t iflash=0; iflash<flashList.second.size(); iflash++) {
328 auto const& flash = flashList.second[iflash];
329 if(flash->TotalPE()<fFlashPeThresh){
333 double tflash = flash->Time();
335 TVector3 rflash(0,flash->YCenter(),flash->ZCenter());
336 TVector3 vdiff = rcrt-rflash;
338 if(
abs(tcrt-tflash)<
abs(tdiff)) {
339 peflash = flash->TotalPE();
343 xyzt.push_back(rflash.X());
344 xyzt.push_back(rflash.Y());
345 xyzt.push_back(rflash.Z());
346 xyzt.push_back(tflash);
348 matchtpc = flashList.first;
350 vector<art::Ptr<recob::OpHit>> hits = findManyHits.at(iflash);
351 for(
auto const&
hit : hits) {
352 double tPmt =
hit->PeakTime();
353 if( tPmt < flashHitT) {
355 flashHitPE =
hit->PE();
358 fGeometryService->OpDetGeoFromOpChannel(
hit->OpChannel()).GetCenter(pos);
359 flashHitxyzt.clear();
360 for(
int i=0; i<3; i++) flashHitxyzt.push_back(pos[i]);
361 flashHitxyzt.push_back(flashHitT);
364 TVector3 rflashHit(pos[0],pos[1],pos[2]);
365 TVector3 vdiffHit = rcrt-rflashHit;
366 flashHitDiff = vdiffHit.Mag();
374 for(
int i=0; i<4; i++) xyzt.push_back(DBL_MAX);
377 fMatchFlash.push_back(matched);
378 fTofFlash.push_back(tdiff);
379 fTofPeFlash.push_back(peflash);
380 fTofXYZTFlash.push_back(xyzt);
381 fDistFlash.push_back(rdiff);
382 fTofTpcFlash.push_back(matchtpc);
383 fTofFlashHit.push_back(tcrt-flashHitT);
384 fTofPeFlashHit.push_back(flashHitPE);
385 fTofXYZTFlashHit.push_back(flashHitxyzt);
386 fDistFlashHit.push_back(flashHitDiff);
396 for(
auto const&
hit : opHitList) {
397 double thit =
hit->PeakTime();
398 if(
hit->PE()<fHitPeThresh){
404 if(
abs(tcrt-thit)<fCoinWindow &&
hit->PE()>pemax) {
409 fGeometryService->OpDetGeoFromOpChannel(
hit->OpChannel()).GetCenter(pos);
412 TVector3 rhit (pos[0],pos[1],pos[2]);
413 TVector3 vdiff = rcrt-rhit;
418 for(
int i=0; i<3; i++) xyzt.push_back(pos[i]);
419 xyzt.push_back(thit);
422 std::cout <<
"thit: "<<thit <<
" , tcrt: " << tcrt <<
" , tdiff " << tdiff << std::endl;
430 for(
int i=0; i<4; i++) xyzt.push_back(DBL_MAX);
432 fMatchHit.push_back(matched);
433 fTofHit.push_back(tdiff);
434 fTofPeHit.push_back(peflash);
435 fTofXYZTHit.push_back(xyzt);
436 fDistHit.push_back(rdiff);
437 fTofTpcHit.push_back(matchtpc);
452 fTofFlashHit.clear();
455 fDistFlashHit.clear();
458 fTofPeFlashHit.clear();
460 fTofXYZTFlash.clear();
461 fTofXYZTFlashHit.clear();
463 fTofTpcFlash.clear();
470 if(hit1->ts1_ns != hit2->ts1_ns)
return false;
471 if(hit1->plane != hit2->plane)
return false;
472 if(hit1->x_pos != hit2->x_pos)
return false;
473 if(hit1->y_pos != hit2->y_pos)
return false;
474 if(hit1->z_pos != hit2->z_pos)
return false;
475 if(hit1->x_err != hit2->x_err)
return false;
476 if(hit1->y_err != hit2->y_err)
return false;
477 if(hit1->z_err != hit2->z_err)
return false;
478 if(hit1->tagger != hit2->tagger)
return false;
vector< vector< double > > fTofXYZTFlashHit
uint64_t m_gate_start_timestamp
vector< int > fTofTpcFlash
art::InputTag fOpHitModuleLabel
process_name opflashCryo1 flashfilter analyze
vector< double > fTofFlash
int fEvent
number of the event being processed
Utilities related to art service access.
vector< bool > fTrackFilt
vector< double > fDistFlash
vector< double > fDistFlashHit
art::InputTag fOpFlashModuleLabel1
vector< bool > fMatchFlash
vector< double > fTofPeFlashHit
map< int, art::InputTag > fFlashLabels
vector< vector< double > > fTofXYZTHit
uint64_t m_trigger_timestamp
vector< double > fTofPeHit
std::string bitName(triggerSource bit)
Returns a mnemonic short name of the beam type.
vector< double > fDistHit
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Access the description of detector geometry.
vector< vector< double > > fTofXYZTFlash
vector< double > fTofPeFlash
Description of geometry of one entire detector.
bool HitCompare(const art::Ptr< CRTHit > &h1, const art::Ptr< CRTHit > &h2)
vector< vector< double > > fCrtXYZErr
int fSubRun
number of the sub-run being processed
CRTPMTMatchingAna(fhicl::ParameterSet const &p)
vector< vector< double > > fCrtXYZT
int fRun
number of the run being processed
vector< double > fTofFlashHit
art::InputTag fTriggerLabel
art::InputTag fCrtHitModuleLabel
art::InputTag fOpFlashModuleLabel2
static double tdiff(const art::Timestamp &ts1, const art::Timestamp &ts2)
triggerSource
Type of beam or beam gate or other trigger source.
art::InputTag fOpFlashModuleLabel3
art::ServiceHandle< art::TFileService > tfs
geo::GeometryCore const * fGeometryService
pointer to Geometry provider
uint64_t m_trigger_gate_diff
void analyze(art::Event const &e) override
art::InputTag fOpFlashModuleLabel0
art framework interface to geometry description
BEGIN_PROLOG could also be cout