315 std::vector<art::Ptr<recob::Track> > tracklist;
316 auto trackListHandle = evt.
getHandle< std::vector<recob::Track> >(
"pandoraTrack");
317 if (trackListHandle){
322 std::vector<art::Ptr<recob::PFParticle> > pfplist;
324 auto PFPListHandle = evt.
getHandle< std::vector<recob::PFParticle> >(
"pandora");
327 std::vector<art::Ptr<recob::Hit>> hitlist;
331 art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, evt,
fTrackModuleLabel);
332 art::FindManyP<anab::T0> trk_t0_assn_v(PFPListHandle, evt ,
"pandora");
333 art::FindManyP<recob::PFParticle> pfp_trk_assn(trackListHandle,evt,
"pandoraTrack");
334 art::FindManyP<anab::T0> fmT0(trackListHandle, evt ,
"pmtrack");
348 std::vector<const sim::SimChannel*> fSimChannels;
350 evt.
getView(
"largeant", fSimChannels);
355 art::FindManyP<anab::T0> fmt0crt2(trackListHandle, evt,
"crtreco");
356 art::FindManyP<anab::CosmicTag> fmctcrt2(trackListHandle, evt,
"crtreco");
373 for(
int k=0;
k < 6;
k++)
377 for(
int k=0;
k < 6;
k++)
383 std::vector< float> peakT_0; std::vector< float> peakT_1; std::vector< float> peakT_2;
384 std::vector< int > tpc_0; std::vector< int > tpc_1; std::vector<int> tpc_2;
385 std::vector< int > wire_0; std::vector< int > wire_1; std::vector<int> wire_2;
386 std::vector< float > int_0; std::vector< float > int_1; std::vector<float> int_2;
387 std::vector< float > amp_0; std::vector< float > amp_1; std::vector<float> amp_2;
388 std::vector<float> rms_0; std::vector<float> rms_1; std::vector<float> rms_2;
389 std::vector<float> hitx_0; std::vector<float> hitx_1; std::vector<float> hitx_2;
390 std::vector<float> hity_0; std::vector<float> hity_1; std::vector<float> hity_2;
391 std::vector<float> hitz_0; std::vector<float> hitz_1; std::vector<float> hitz_2;
392 std::vector<float> hitz_wire2; std::vector<float> startcosxyz; std::vector<float> endcosxyz;
393 std::vector<float> dT_buffer; std::vector<float> gof, multi,truerms;
398 size_t NTracks = tracklist.size();
399 for(
size_t i=0; i<NTracks;++i){
403 peakT_0.clear(); peakT_1.clear(); peakT_2.clear();
404 tpc_0.clear(); tpc_1.clear(); tpc_2.clear();
405 wire_0.clear(); wire_1.clear(); wire_2.clear();
406 int_0.clear(); int_1.clear() ; int_2.clear();
407 amp_0.clear(); amp_1.clear() ; amp_2.clear();
408 rms_0.clear(); rms_1.clear(); rms_2.clear();
409 hitx_0.clear(); hitx_1.clear(); hitx_2.clear();
410 hity_0.clear(); hity_1.clear(); hity_2.clear();
411 hitz_0.clear(); hitz_1.clear(); hitz_2.clear();
412 hitz_wire2.clear(); startcosxyz.clear();endcosxyz.clear();
413 dT_buffer.clear(); gof.clear(); multi.clear();
420 double t_zero=-999999;
447 double this_t0crt2=-DBL_MAX;
449 double this_crt2x0 = -DBL_MAX;
450 double this_crt2x1 = -DBL_MAX;
451 double this_crt2y0 = -DBL_MAX;
452 double this_crt2y1 = -DBL_MAX;
453 double this_crt2z0 = -DBL_MAX;
454 double this_crt2z1 = -DBL_MAX;
458 if(fmt0crt2.isValid()){
459 auto const& vt0crt2 = fmt0crt2.at(i);
460 if (!vt0crt2.empty()){
461 this_t0crt2 = vt0crt2[0]->Time();
467 if (fmctcrt2.isValid()){
468 auto const& vctcrt2 = fmctcrt2.at(i);
469 if (!vctcrt2.empty()){
470 this_crt2x0 = vctcrt2[0]->EndPoint1()[0];
471 this_crt2x1 = vctcrt2[0]->EndPoint2()[0];
472 this_crt2y0 = vctcrt2[0]->EndPoint1()[1];
473 this_crt2y1 = vctcrt2[0]->EndPoint2()[1];
474 this_crt2z0 = vctcrt2[0]->EndPoint1()[2];
475 this_crt2z1 = vctcrt2[0]->EndPoint2()[2];
499 double theta_xz = std::atan2(dir_start.X(), dir_start.Z());
500 double theta_yz = std::atan2(dir_start.Y(), dir_start.Z());
506 auto allHits=fmthm.at(i);
507 double ticksoffset=0;
510 if (this_t0crt2 > -DBL_MAX) ticksoffset = this_t0crt2/500.+detProp.GetXTicksOffset (allHits[0]->
WireID().
Plane, allHits[0]->
WireID().
TPC, allHits[0]->
WireID().Cryostat);
512 xoffset = detProp.ConvertTicksToX(ticksoffset,allHits[0]->
WireID());
513 std::cout<<
"tickoffset , x offset "<<ticksoffset<<
" "<<xoffset<<
" default term "<<detProp.GetXTicksOffset (allHits[10]->
WireID().
Plane, allHits[10]->
WireID().
TPC, allHits[10]->
WireID().Cryostat)<<
std::endl;
519 auto vhit=fmthm.at(i);
520 auto vmeta=fmthm.data(i);
521 for (
size_t ii = 0; ii<vhit.size(); ++ii){
522 bool fBadhit =
false;
528 if (vmeta[ii]->
Index()>=tracklist[i]->NumberTrajectoryPoints()){
529 throw cet::exception(
"Calorimetry_module.cc") <<
"Requested track trajectory index "<<vmeta[ii]->Index()<<
" exceeds the total number of trajectory points "<<tracklist[i]->NumberTrajectoryPoints()<<
" for track index "<<i<<
". Something is wrong with the track reconstruction. Please contact tjyang@fnal.gov!!";
531 if (!tracklist[i]->HasValidPoint(vmeta[ii]->
Index())){
537 auto loc = tracklist[i]->LocationAtPoint(vmeta[ii]->
Index());
543 if (fBadhit)
continue;
544 if (zpos<-100)
continue;
545 planenum=vhit[ii]->WireID().Plane;
547 peakT_0.push_back(vhit[ii]->PeakTime());
550 int_0.push_back(vhit[ii]->Integral());
551 amp_0.push_back(vhit[ii]->PeakAmplitude());
552 rms_0.push_back(vhit[ii]->RMS());
553 hitx_0.push_back(xpos);
554 hity_0.push_back(ypos);
555 hitz_0.push_back(zpos);
558 peakT_1.push_back(vhit[ii]->PeakTime());
561 int_1.push_back(vhit[ii]->Integral());
562 amp_1.push_back(vhit[ii]->PeakAmplitude());
563 rms_1.push_back(vhit[ii]->RMS());
564 hitx_1.push_back(xpos);
565 hity_1.push_back(ypos);
566 hitz_1.push_back(zpos);
570 bool removehit=
false;
571 for (
size_t i2=0; i2<hitlist.size(); ++i2) {
572 if (vhit[ii].
key() == hitlist[i2].key())
continue;
573 if (vhit[ii]->Channel()!=hitlist[i2]->Channel())
continue;
574 if ((vhit[ii]->PeakTime()+30<hitlist[i2]->PeakTime()-30) || (vhit[ii]->PeakTime()-30>hitlist[i2]->PeakTime()+30))
continue;
578 if (removehit)
continue;
581 peakT_2.push_back(vhit[ii]->PeakTime()-this_t0crt2/500.0);
584 int_2.push_back(vhit[ii]->Integral());
585 amp_2.push_back(vhit[ii]->PeakAmplitude());
586 rms_2.push_back(vhit[ii]->RMS());
587 dT_buffer.push_back(vhit[ii]->EndTick()-vhit[ii]->StartTick());
588 hitx_2.push_back(xpos);
589 hity_2.push_back(ypos);
590 hitz_2.push_back(zpos);
591 gof.push_back(vhit[ii]->GoodnessOfFit());
592 multi.push_back(vhit[ii]->Multiplicity());
595 unsigned int wireno=vhit[ii]->WireID().Wire;
597 hitz_wire2.push_back(xyzStart[2]);
601 unsigned int t0=vhit[ii]->PeakTime()-3*(vhit[ii]->RMS());
603 unsigned int t1=vhit[ii]->PeakTime()+3*(vhit[ii]->RMS());
604 if(t1>6000) t1=6000-1;
716 for(
size_t sc=0;sc<fSimChannels.size();++sc){
717 if(fSimChannels[sc]->Channel()==vhit[ii]->Channel()) chan=fSimChannels[sc];
720 const auto &tdcidemap = chan->
TDCIDEMap();
724 double hit_rmsvalue=0;
725 for(
auto mapitr = tdcidemap.begin(); mapitr != tdcidemap.end(); mapitr++){
726 if(!((*mapitr).first>t0 && (*mapitr).first<=t1))
continue;
729 const std::vector<sim::IDE> idevec = (*mapitr).second;
730 for(
size_t iv = 0; iv < idevec.size(); ++iv){
731 hit_nelec+=idevec[iv].numElectrons;
735 int j=(*mapitr).first;
736 mean_t+=double(j)*double(hit_nelec);
737 double jterm=double(j)*double(j)*double(hit_nelec);
738 mean_t2=mean_t2+jterm;
742 hit_rmsvalue=sqrt(mean_t2-mean_t*mean_t);
744 truermsb=hit_rmsvalue;
747 truerms.push_back(truermsb);
754 if(peakT_2.size()<10)
continue;
755 max_value=*std::max_element(peakT_2.begin(),peakT_2.end());
756 min_value=*std::min_element(peakT_2.begin(),peakT_2.end());
758 std::cout<<max_value<<
" "<<min_value<<
std::endl;
806 startcosxyz.push_back(dir_start.X());
807 startcosxyz.push_back(dir_start.Y());
808 startcosxyz.push_back(dir_start.Z());
809 endcosxyz.push_back(dir_end.X());
810 endcosxyz.push_back(dir_end.Y());
811 endcosxyz.push_back(dir_end.Z());
830 t0crt2.push_back(this_t0crt2/500.0);
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
code to link reconstructed objects back to the MC truth information
std::vector< std::vector< float > > goodnessoffit2
std::vector< std::vector< float > > trkdq_amp0
constexpr std::uint32_t timeLow() const
std::vector< float > trackthetayz
std::vector< float > crtreco_y0
std::vector< float > trklen
Energy deposited on a readout channel by simulated tracks.
std::vector< float > crtreco_z1
Handle< PROD > getHandle(SelectorBase const &) const
std::vector< std::vector< float > > hit_rms2
std::vector< std::vector< float > > trkhitx1
std::vector< std::vector< int > > hit_wire2
std::vector< float > trkendy
constexpr std::uint32_t timeHigh() const
std::vector< float > trkstartx
std::vector< int > tot_trks
std::vector< std::vector< float > > trkhity0
Vector_t VertexDirection() const
std::vector< std::vector< float > > trkdq_amp2
art::InputTag fRawProducerLabel
std::vector< std::vector< float > > hit_rms0
std::vector< std::vector< int > > hit_tpc1
std::vector< float > trackthetaxz
std::vector< std::vector< float > > trkhitz1
std::vector< double > T0_values
std::vector< std::vector< float > > hit_deltaT
geo::GeometryCore const * fGeometry
std::vector< std::vector< float > > trkdq_int2
std::vector< double > crt2tickoffset
std::vector< float > crtreco_z0
std::vector< std::vector< float > > multiplicity2
std::vector< float > trkstarty
std::vector< std::vector< int > > hit_tpc0
std::vector< float > trkendx_crt2
double Length(size_t p=0) const
Access to various track properties.
std::vector< std::vector< float > > hit_rms1
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
Point_t const & Vertex() const
std::string fTrackModuleLabel
std::vector< std::vector< float > > hit_peakT0
std::vector< float > trkendx
std::vector< std::vector< int > > hit_tpc2
std::vector< std::vector< float > > trkhitx2
SubRunNumber_t subRun() const
std::vector< double > t0crt2
static int max(int a, int b)
std::vector< std::vector< float > > trkhitx0
std::vector< std::vector< int > > hit_wire1
std::vector< float > xprojectedlen
std::vector< float > trkendz
std::vector< std::vector< float > > trkhitz0
std::vector< std::vector< float > > trkhitz_wire2
std::vector< std::vector< float > > trkhity1
std::vector< std::vector< float > > trkdq_int1
std::vector< std::vector< float > > trkhitz2
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
std::vector< std::vector< float > > trkendcosxyz
std::vector< float > trkstartz
int GetNActiveFembsForAPA(art::Event const &evt, int apa) const
Get number of active fembs in an APA.
std::string fHitsModuleLabel
std::vector< float > crtreco_x0
Vector_t EndDirection() const
std::vector< std::vector< float > > hit_peakT1
std::vector< std::vector< int > > hit_wire0
std::vector< std::vector< float > > trkdq_int0
void WireEndPoints(geo::WireID const &wireid, double *xyzStart, double *xyzEnd) const
Fills two arrays with the coordinates of the wire end points.
EventNumber_t event() const
Point_t const & End() const
std::vector< std::vector< float > > hit_rms_true
std::vector< float > crtreco_x1
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
std::vector< std::vector< float > > trkhity2
recob::tracking::Plane Plane
std::vector< float > crtreco_y1
std::vector< std::vector< float > > trkstartcosxyz
ProtoDUNEDataUtils fDataUtils
std::vector< float > trkstartx_crt2
std::vector< std::vector< float > > trkdq_amp1
std::vector< std::vector< float > > hit_peakT2
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
cet::coded_exception< error, detail::translate > exception
QTextStream & endl(QTextStream &s)