221 std::cout<<
"============================================"<<std::endl
222 <<
"Run = "<<
event.run()<<
", SubRun = "<<
event.subRun()<<
", Event = "<<
event.id().event()<<std::endl
223 <<
"============================================"<<std::endl;
230 auto particleHandle =
event.getValidHandle<std::vector<simb::MCParticle>>(
fSimModuleLabel);
233 auto crtHitHandle =
event.getValidHandle<std::vector<sbn::crt::CRTHit>>(
fCRTHitLabel);
236 auto crtTrackHandle =
event.getValidHandle<std::vector<sbn::crt::CRTTrack>>(
fCRTTrackLabel);
239 art::FindManyP<sbn::crt::CRTHit> findManyHits(crtTrackHandle, event,
fCRTTrackLabel);
244 std::map<int, std::vector<sbn::crt::CRTTrack>> crtTracks;
247 for(
auto const&
track : (*crtTrackHandle)){
250 if(trueId == -99999)
continue;
251 crtTracks[trueId].push_back(
track);
254 std::map<int, std::vector<sbn::crt::CRTHit>> crtHits;
255 double minHitTime = 99999;
256 double maxHitTime = -99999;
258 for(
auto const&
hit : (*crtHitHandle)){
259 double hitTime = (double)(
int)
hit.ts1_ns * 1
e-3;
260 if(hitTime < minHitTime) minHitTime = hitTime;
261 if(hitTime > maxHitTime) maxHitTime = hitTime;
265 if(trueId == -99999)
continue;
266 crtHits[trueId].push_back(
hit);
273 std::map<int, simb::MCParticle> particles;
274 for (
auto const& particle: (*particleHandle)){
275 int partId = particle.TrackId();
276 particles[partId] = particle;
279 double time = particle.T() * 1
e-3;
280 if(time < minHitTime || time > maxHitTime)
continue;
282 if(!(
std::abs(particle.PdgCode()) == 13 && particle.Mother()==0))
continue;
285 if(crtHits.find(partId) == crtHits.end())
continue;
287 std::vector<std::string> hitTaggers;
288 for(
auto const&
hit : crtHits[partId]){
289 if(std::find(hitTaggers.begin(), hitTaggers.end(),
hit.tagger) != hitTaggers.end())
continue;
290 hitTaggers.push_back(
hit.tagger);
293 if(nPlanesHit < 2)
continue;
295 double momentum = particle.P();
296 TVector3 start (particle.Vx(), particle.Vy(), particle.Vz());
297 TVector3
end (particle.EndX(), particle.EndY(), particle.EndZ());
298 double theta = (
end-start).Theta();
299 double phi = (
end-start).Phi();
300 for(
auto const&
tagger : hitTaggers){
309 if(crtTracks.find(partId) == crtTracks.end())
continue;
310 for(
auto const&
tagger : hitTaggers){
324 for(
auto const&
track : (*crtTrackHandle)){
334 std::vector<art::Ptr<sbn::crt::CRTHit>> hits = findManyHits.at(track_i);
338 if(particles.find(trueId) == particles.end())
continue;
343 double denominator = (
end - start).Mag();
346 simb::MCParticle particle = particles[trueId];
348 int nTraj = particle.NumberTrajectoryPoints();
351 for(
int j = 0; j < nTraj; j++){
352 TVector3 pt (particle.Vx(j), particle.Vy(j), particle.Vz(j));
353 if (pt.X() >= crtLims[0] && pt.X() <= crtLims[3] && pt.Y() >= crtLims[1] && pt.Y() <= crtLims[4] && pt.Z() >= crtLims[2] && pt.Z() <= crtLims[5]){
354 double numerator = ((pt - start).Cross(pt -
end)).Mag();
355 aveDCA += numerator/denominator;
360 aveDCA = aveDCA/npts;
364 for(
auto const&
hit : hits){
366 if(trueCross.X() == -99999)
continue;
368 double dist = std::sqrt(std::pow(
hit->x_pos - trueCross.X(), 2)
369 + std::pow(
hit->y_pos - trueCross.Y(), 2)
370 + std::pow(
hit->z_pos - trueCross.Z(), 2));
386 auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
int TrueIdFromTrackId(const art::Event &event, int track_i)
void Draw(detinfo::DetectorClocksData const &clockData, const art::Event &event)
art::InputTag fCRTHitLabel
name of CRT hit producer
std::map< std::string, TH1D * > hCrossDistance
art::InputTag fSimModuleLabel
name of detsim producer
process_name use argoneut_mc_hitfinder track
std::map< std::string, TH1D * > hEffPhiReco
int TrueIdFromHitId(const art::Event &event, int hit_i)
void Initialize(const art::Event &event)
CRTBackTracker fCrtBackTrack
art::InputTag fCRTTrackLabel
name of CRT track producer
int fPlotTrackID
id of track to plot
std::map< std::string, TH1D * > hTrackDist
auto end(FixedBins< T, C > const &) noexcept
std::map< std::string, TH1D * > hEffThetaReco
std::map< std::string, TH1D * > hEffMomTotal
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
geo::Point_t TaggerCrossingPoint(std::string taggerName, const simb::MCParticle &particle)
void SetDrawTrueTracks(bool tf)
bool fVeryVerbose
print more information about what's going on
std::map< std::string, TH1D * > hEffPhiTotal
std::vector< double > CRTLimits() const
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
bool fVerbose
print information about what's going on
std::map< std::string, TH1D * > hEffThetaTotal
std::map< std::string, TH1D * > hEffMomReco
BEGIN_PROLOG could also be cout