225 std::cout<<
"============================================"<<std::endl
226 <<
"Run = "<<
event.run()<<
", SubRun = "<<
event.subRun()<<
", Event = "<<
event.id().event()<<std::endl
227 <<
"============================================"<<std::endl;
235 art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
236 auto particleHandle =
event.getValidHandle<std::vector<simb::MCParticle>>(
fSimModuleLabel);
239 art::Handle< std::vector<sbn::crt::CRTHit>> crtHitHandle;
240 std::vector<art::Ptr<sbn::crt::CRTHit> > crtHitList;
242 art::fill_ptr_vector(crtHitList, crtHitHandle);
245 auto tpcTrackHandle =
event.getValidHandle<std::vector<recob::Track>>(
fTPCTrackLabel);
246 art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event,
fTPCTrackLabel);
254 std::map<int, simb::MCParticle> particles;
256 for (
auto const& particle: (*particleHandle)){
259 int partID = particle.TrackId();
260 particles[partID] = particle;
264 std::map<int, std::vector<std::string>> crtTaggerMap;
265 std::vector<sbn::crt::CRTHit> crtHits;
267 double minHitTime = 99999;
268 double maxHitTime = -99999;
270 for(
auto const&
hit : (*crtHitHandle)){
271 double hitTime = (double)(
int)
hit.ts1_ns * 1
e-3;
272 if(hitTime < minHitTime) minHitTime = hitTime;
273 if(hitTime > maxHitTime) maxHitTime = hitTime;
275 crtHits.push_back(
hit);
278 if(trueId == -99999)
continue;
279 if(crtTaggerMap.find(trueId) == crtTaggerMap.end()){
280 crtTaggerMap[trueId].push_back(
hit.tagger);
282 else if(std::find(crtTaggerMap[trueId].
begin(), crtTaggerMap[trueId].
end(),
hit.tagger) == crtTaggerMap[trueId].end()){
283 crtTaggerMap[trueId].push_back(
hit.tagger);
292 auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
293 auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(event, clockData);
296 for (
auto const& tpcTrack : (*tpcTrackHandle)){
298 std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
300 if(particles.find(trackTrueID) == particles.end())
continue;
302 if(!(
std::abs(particles[trackTrueID].PdgCode()) == 13 && particles[trackTrueID].Mother() == 0))
continue;
305 double trueTime = particles[trackTrueID].T() * 1
e-3;
306 if(trueTime < minHitTime || trueTime > maxHitTime)
continue;
312 double sin_angle = -99999;
313 if(closest.dca != -99999){
314 hDCA[closest.thishit.tagger]->Fill(closest.dca);
315 hDCA[
"All"]->Fill(closest.dca);
316 sin_angle = closest.dca/closest.extrapLen;
317 hDoL[closest.thishit.tagger]->Fill(sin_angle);
318 hDoL[
"All"]->Fill(sin_angle);
319 hT0[closest.thishit.tagger]->Fill(closest.t0);
320 hT0[
"All"]->Fill(closest.t0);
325 if(hitTrueID == trackTrueID && hitTrueID != -99999){
326 hMatchDCA[closest.thishit.tagger]->Fill(closest.dca);
328 hMatchDoL[closest.thishit.tagger]->Fill(sin_angle);
330 hMatchT0[closest.thishit.tagger]->Fill(closest.t0);
334 hNoMatchDCA[closest.thishit.tagger]->Fill(closest.dca);
336 hNoMatchDoL[closest.thishit.tagger]->Fill(sin_angle);
338 hNoMatchT0[closest.thishit.tagger]->Fill(closest.t0);
344 for(
int i = 0; i < nbins; i++){
345 double DCAcut =
hEffDCATotal.begin()->second->GetBinCenter(i);
348 if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
349 for(
auto const&
tagger : crtTaggerMap[trackTrueID]){
353 if(closest.dca < DCAcut && closest.dca != -99999){
359 if(closest.dca < DCAcut && closest.dca != -99999){
365 if(closest.dca < DCAcut && closest.dca != -99999){
370 double hitTime = closest.thishit.ts1_ns * 1
e-3;
371 if(particles.find(trackTrueID) != particles.end()){
372 if(
std::abs(hitTime - trueTime) < 2.){
383 for(
int i = 0; i < nbins; i++){
384 double DCAcut =
hEffDoLTotal.begin()->second->GetBinCenter(i);
387 if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
388 for(
auto const&
tagger : crtTaggerMap[trackTrueID]){
392 if(sin_angle < DCAcut && closest.dca != -99999){
398 if(sin_angle < DCAcut && closest.dca != -99999){
404 if(sin_angle < DCAcut && closest.dca != -99999){
409 double hitTime = closest.thishit.ts1_ns * 1
e-3;
410 if(particles.find(trackTrueID) != particles.end()){
411 if(
std::abs(hitTime - trueTime) < 2.){
420 double fixedCut = 30.;
423 if(crtTaggerMap.find(trackTrueID) != crtTaggerMap.end()){
424 for(
auto const&
tagger : crtTaggerMap[trackTrueID]){
428 if(closest.dca < fixedCut && closest.dca >=0 ){
434 if(closest.dca < fixedCut && closest.dca >=0){
440 if(closest.dca < fixedCut && closest.dca >= 0){
445 double hitTime = closest.thishit.ts1_ns * 1
e-3;
446 if(particles.find(trackTrueID) != particles.end()){
447 double trueTime = particles[trackTrueID].T() * 1
e-3;
448 if(
std::abs(hitTime - trueTime) < 2.){
art::InputTag fSimModuleLabel
name of detsim producer
std::map< std::string, TH1D * > hMatchT0
std::map< std::string, TH1D * > hPurityLengthTotal
std::map< std::string, TH1D * > hEffDCAReco
matchCand GetClosestCRTHit(detinfo::DetectorPropertiesData const &detProp, recob::Track tpcTrack, std::pair< double, double > t0MinMax, std::vector< sbn::crt::CRTHit > crtHits, int driftDirection)
std::map< std::string, TH1D * > hEffDoLReco
std::map< std::string, TH1D * > hDoL
art::InputTag fCRTHitLabel
name of CRT producer
std::map< std::string, TH1D * > hPurityLengthReco
std::map< std::string, TH1D * > hEffLengthReco
std::map< std::string, TH1D * > hPurityDCAReco
std::map< std::string, TH1D * > hNoMatchDoL
std::map< std::string, TH1D * > hDCA
int TrueIdFromTotalEnergy(const art::Event &event, const sbnd::crt::CRTData &data)
int TrueIdFromHitId(const art::Event &event, int hit_i)
std::map< std::string, TH1D * > hPurityDoLReco
void Initialize(const art::Event &event)
std::map< std::string, TH1D * > hNoMatchT0
std::map< std::string, TH1D * > hMatchDoL
std::map< std::string, TH1D * > hPurityDCATotal
auto end(FixedBins< T, C > const &) noexcept
std::map< std::string, TH1D * > hT0
int TrueParticleIDFromTotalRecoHits(detinfo::DetectorClocksData const &clockData, const std::vector< art::Ptr< recob::Hit > > &hits, bool rollup_unsaved_ids=1)
art::InputTag fTPCTrackLabel
name of CRT producer
std::map< std::string, TH1D * > hEffLengthTotal
std::map< std::string, TH1D * > hMatchDCA
std::map< std::string, TH1D * > hEffDoLTotal
auto begin(FixedBins< T, C > const &) noexcept
std::map< std::string, TH1D * > hPurityDoLTotal
bool fVerbose
print information about what's going on
std::map< std::string, TH1D * > hNoMatchDCA
CRTBackTracker fCrtBackTrack
BEGIN_PROLOG could also be cout
std::map< std::string, TH1D * > hEffDCATotal