26 #include "nusimdata/SimulationBase/MCParticle.h"
30 #include "art/Framework/Core/EDAnalyzer.h"
31 #include "art/Framework/Principal/Event.h"
32 #include "art/Framework/Principal/Handle.h"
33 #include "art/Framework/Services/Registry/ServiceHandle.h"
34 #include "art_root_io/TFileService.h"
35 #include "art/Framework/Core/ModuleMacros.h"
36 #include "canvas/Persistency/Common/FindManyP.h"
37 #include "canvas/Utilities/Exception.h"
42 #include "messagefacility/MessageLogger/MessageLogger.h"
43 #include "fhiclcpp/ParameterSet.h"
44 #include "fhiclcpp/types/Table.h"
45 #include "fhiclcpp/types/Atom.h"
46 #include "cetlib/pow.h"
72 Name(
"SimModuleLabel"),
73 Comment(
"tag of detector simulation data product")
78 Comment(
"tag of CRT hit producer data product")
82 Name(
"TPCTrackLabel"),
83 Comment(
"tag of tpc track producer data product")
88 Comment(
"Print information about what's going on")
110 virtual void analyze(
const art::Event& event)
override;
113 virtual void endJob()
override;
134 std::map<std::string, TH1D*>
hDCA;
138 std::map<std::string, TH1D*>
hDoL;
142 std::map<std::string, TH1D*>
hT0;
166 , fSimModuleLabel (config().SimModuleLabel())
167 , fCRTHitLabel (config().CRTHitLabel())
169 , fVerbose (config().
Verbose())
170 , t0Alg (config().CRTT0Alg())
171 , fCrtBackTrack (config().CrtBackTrack())
181 art::ServiceHandle<art::TFileService>
tfs;
183 std::string
tagger =
"All";
187 hDCA[tagger] = tfs->make<TH1D>(Form(
"DCA_%s", tagger.c_str()),
"", 50, 0, 100);
188 hMatchDCA[tagger] = tfs->make<TH1D>(Form(
"MatchDCA_%s", tagger.c_str()),
"", 50, 0, 100);
189 hNoMatchDCA[tagger] = tfs->make<TH1D>(Form(
"NoMatchDCA_%s", tagger.c_str()),
"", 50, 0, 100);
191 hDoL[tagger] = tfs->make<TH1D>(Form(
"DoL_%s", tagger.c_str()),
"", 100, 0, 0.25);
192 hMatchDoL[tagger] = tfs->make<TH1D>(Form(
"MatchDoL_%s", tagger.c_str()),
"", 100, 0, 0.25);
193 hNoMatchDoL[tagger] = tfs->make<TH1D>(Form(
"NoMatchDoL_%s", tagger.c_str()),
"", 100, 0, 0.25);
195 hT0[tagger] = tfs->make<TH1D>(Form(
"T0_%s", tagger.c_str()),
"", 600, -3000, 3000);
196 hMatchT0[tagger] = tfs->make<TH1D>(Form(
"MatchT0_%s", tagger.c_str()),
"", 600, -3000, 3000);
197 hNoMatchT0[tagger] = tfs->make<TH1D>(Form(
"NoMatchT0_%s", tagger.c_str()),
"", 600, -3000, 3000);
199 hEffDCATotal[tagger] = tfs->make<TH1D>(Form(
"EffDCATotal_%s", tagger.c_str()),
"", 50, 0, 100);
200 hEffDCAReco[tagger] = tfs->make<TH1D>(Form(
"EffDCAReco_%s", tagger.c_str()),
"", 50, 0, 100);
201 hEffDoLTotal[tagger] = tfs->make<TH1D>(Form(
"EffDoLTotal_%s", tagger.c_str()),
"", 100, 0, 0.25);
202 hEffDoLReco[tagger] = tfs->make<TH1D>(Form(
"EffDoLReco_%s", tagger.c_str()),
"", 100, 0, 0.25);
203 hEffLengthTotal[tagger] = tfs->make<TH1D>(Form(
"EffLengthTotal_%s", tagger.c_str()),
"", 20, 0, 600);
204 hEffLengthReco[tagger] = tfs->make<TH1D>(Form(
"EffLengthReco_%s", tagger.c_str()),
"", 20, 0, 600);
206 hPurityDCATotal[tagger] = tfs->make<TH1D>(Form(
"PurityDCATotal_%s", tagger.c_str()),
"", 50, 0, 100);
207 hPurityDCAReco[tagger] = tfs->make<TH1D>(Form(
"PurityDCAReco_%s", tagger.c_str()),
"", 50, 0, 100);
208 hPurityDoLTotal[tagger] = tfs->make<TH1D>(Form(
"PurityDoLTotal_%s", tagger.c_str()),
"", 100, 0, 0.25);
209 hPurityDoLReco[tagger] = tfs->make<TH1D>(Form(
"PurityDoLReco_%s", tagger.c_str()),
"", 100, 0, 0.25);
210 hPurityLengthTotal[tagger] = tfs->make<TH1D>(Form(
"PurityLengthTotal_%s", tagger.c_str()),
"", 20, 0, 600);
211 hPurityLengthReco[tagger] = tfs->make<TH1D>(Form(
"PurityLengthReco_%s", tagger.c_str()),
"", 20, 0, 600);
215 if(
fVerbose)
std::cout<<
"----------------- CRT T0 Matching Ana Module -------------------"<<std::endl;
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){
315 hDCA[
"All"]->Fill(closest.
dca);
318 hDoL[
"All"]->Fill(sin_angle);
320 hT0[
"All"]->Fill(closest.
t0);
325 if(hitTrueID == trackTrueID && hitTrueID != -99999){
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){
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){
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){
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
Functions to help with numbers.
std::map< std::string, TH1D * > hDoL
art::InputTag fCRTHitLabel
name of CRT producer
double ts1_ns
Timestamp T1 ([signal time w.r.t. Trigger time]), in UTC absolute time scale in nanoseconds from the ...
CRTTaggerGeo GetTagger(std::string taggerName) const
std::map< std::string, TH1D * > hPurityLengthReco
std::map< std::string, TH1D * > hEffLengthReco
Declaration of signal hit object.
fhicl::Atom< bool > Verbose
std::map< std::string, TH1D * > hPurityDCAReco
fhicl::Atom< art::InputTag > TPCTrackLabel
std::map< std::string, TH1D * > hNoMatchDoL
virtual void analyze(const art::Event &event) override
std::map< std::string, TH1D * > hDCA
CRTT0MatchingAna(Parameters const &config)
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
fhicl::Table< CRTBackTracker::Config > CrtBackTrack
void Initialize(const art::Event &event)
fhicl::Atom< art::InputTag > CRTHitLabel
virtual void beginJob() override
std::map< std::string, TH1D * > hNoMatchT0
double DistToCrtHit(TVector3 trackPos, sbn::crt::CRTHit crtHit)
fhicl::Atom< art::InputTag > SimModuleLabel
std::map< std::string, TH1D * > hMatchDoL
virtual void endJob() override
art::EDAnalyzer::Table< Config > Parameters
fhicl::Table< CRTT0MatchAlg::Config > CRTT0Alg
std::map< std::string, TH1D * > hPurityDCATotal
auto end(FixedBins< T, C > const &) noexcept
BEGIN_PROLOG vertical distance to the surface Name
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)
size_t NumTaggers() const
BEGIN_PROLOG don t mess with this TPCTrackLabel
art::InputTag fTPCTrackLabel
name of CRT producer
Definition of data types for geometry description.
std::map< std::string, TH1D * > hEffLengthTotal
Provides recob::Track data product.
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
object containing MC truth information necessary for making RawDigits and doing back tracking ...
stream1 can override from command line with o or output services user sbnd
std::map< std::string, TH1D * > hNoMatchDCA
art::ServiceHandle< art::TFileService > tfs
CRTBackTracker fCrtBackTrack
std::string tagger
Name of the CRT wall (in the form of strings).
BEGIN_PROLOG could also be cout
std::map< std::string, TH1D * > hEffDCATotal