18 #include "art_root_io/TFileService.h" 19 #include "art_root_io/TFileDirectory.h" 21 #include "canvas/Persistency/Common/FindManyP.h" 26 #include "dune/Protodune/singlephase/DataUtils/ProtoDUNEDataUtils.h" 65 #include "TDirectory.h" 73 #include "TGraphErrors.h" 76 #include "TTimeStamp.h" 82 #include "TTimeStamp.h" 86 #include "TLorentzVector.h" 89 #include "TPaveStats.h" 224 fEventTree = tfs->make<TTree>(
"Event",
"Event Tree from Reco");
316 std::vector<art::Ptr<recob::Track> > tracklist;
317 auto trackListHandle = evt.
getHandle< std::vector<recob::Track> >(
"pandoraTrack");
318 if (trackListHandle){
324 std::vector<art::Ptr<recob::PFParticle> > pfplist;
325 auto PFPListHandle = evt.
getHandle< std::vector<recob::PFParticle> >(
"pandora");
328 std::vector<art::Ptr<recob::Hit>> hitlist;
332 art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, evt,
fTrackModuleLabel);
333 art::FindManyP<anab::T0> trk_t0_assn_v(PFPListHandle, evt ,
"pandora");
334 art::FindManyP<recob::PFParticle> pfp_trk_assn(trackListHandle,evt,
"pandoraTrack");
335 art::FindManyP<anab::T0> fmT0(trackListHandle, evt ,
"pmtrack");
338 std::vector<art::Ptr<raw::RawDigit> > rawlist;
343 std::vector<art::Ptr<recob::Wire> > wirelist;
348 std::cout<<
"rawlist size, wirelist size"<<rawlist.size()<<
" "<<wirelist.size()<<
std::endl;
349 std::vector<const sim::SimChannel*> fSimChannels;
351 evt.
getView(
"largeant", fSimChannels);
356 art::FindManyP<anab::T0> fmt0crt2(trackListHandle, evt,
"crtreco");
357 art::FindManyP<anab::CosmicTag> fmctcrt2(trackListHandle, evt,
"crtreco");
374 for(
int k=0;
k < 6;
k++)
378 for(
int k=0;
k < 6;
k++)
384 std::vector< float> peakT_0; std::vector< float> peakT_1; std::vector< float> peakT_2; std::vector<float> hit_peakTrawb;
385 std::vector< int > tpc_0; std::vector< int > tpc_1; std::vector<int> tpc_2;
386 std::vector< int > wire_0; std::vector< int > wire_1; std::vector<int> wire_2;
387 std::vector< float > int_0; std::vector< float > int_1; std::vector<float> int_2;
388 std::vector< float > amp_0; std::vector< float > amp_1; std::vector<float> amp_2;
389 std::vector<float> rms_0; std::vector<float> rms_1; std::vector<float> rms_2, rms_raw2;
390 std::vector<float> hitx_0; std::vector<float> hitx_1; std::vector<float> hitx_2;
391 std::vector<float> hity_0; std::vector<float> hity_1; std::vector<float> hity_2;
392 std::vector<float> hitz_0; std::vector<float> hitz_1; std::vector<float> hitz_2;
393 std::vector<float> hitz_wire2; std::vector<float> startcosxyz; std::vector<float> endcosxyz;
394 std::vector<float> dT_buffer; std::vector<float> gof, multi,truerms;
399 size_t NTracks = tracklist.size();
400 for(
size_t i=0; i<NTracks;++i){
404 peakT_0.clear(); peakT_1.clear(); peakT_2.clear(); hit_peakTrawb.clear();
405 tpc_0.clear(); tpc_1.clear(); tpc_2.clear();
406 wire_0.clear(); wire_1.clear(); wire_2.clear();
407 int_0.clear(); int_1.clear() ; int_2.clear();
408 amp_0.clear(); amp_1.clear() ; amp_2.clear();
409 rms_0.clear(); rms_1.clear(); rms_2.clear();rms_raw2.clear();
410 hitx_0.clear(); hitx_1.clear(); hitx_2.clear();
411 hity_0.clear(); hity_1.clear(); hity_2.clear();
412 hitz_0.clear(); hitz_1.clear(); hitz_2.clear();
413 hitz_wire2.clear(); startcosxyz.clear();endcosxyz.clear();
414 dT_buffer.clear(); gof.clear(); multi.clear();
421 double t_zero=-999999;
448 double this_t0crt2=-DBL_MAX;
450 double this_crt2x0 = -DBL_MAX;
451 double this_crt2x1 = -DBL_MAX;
452 double this_crt2y0 = -DBL_MAX;
453 double this_crt2y1 = -DBL_MAX;
454 double this_crt2z0 = -DBL_MAX;
455 double this_crt2z1 = -DBL_MAX;
459 if(fmt0crt2.isValid()){
460 auto const& vt0crt2 = fmt0crt2.at(i);
461 if (!vt0crt2.empty()){
462 this_t0crt2 = vt0crt2[0]->Time();
468 if (fmctcrt2.isValid()){
469 auto const& vctcrt2 = fmctcrt2.at(i);
470 if (!vctcrt2.empty()){
471 this_crt2x0 = vctcrt2[0]->EndPoint1()[0];
472 this_crt2x1 = vctcrt2[0]->EndPoint2()[0];
473 this_crt2y0 = vctcrt2[0]->EndPoint1()[1];
474 this_crt2y1 = vctcrt2[0]->EndPoint2()[1];
475 this_crt2z0 = vctcrt2[0]->EndPoint1()[2];
476 this_crt2z1 = vctcrt2[0]->EndPoint2()[2];
500 double theta_xz = std::atan2(dir_start.X(), dir_start.Z());
501 double theta_yz = std::atan2(dir_start.Y(), dir_start.Z());
507 auto allHits=fmthm.at(i);
508 double ticksoffset=0;
511 if (this_t0crt2 > -DBL_MAX) ticksoffset = this_t0crt2/500.+detProp.GetXTicksOffset (allHits[0]->
WireID().
Plane, allHits[0]->
WireID().
TPC, allHits[0]->
WireID().Cryostat);
513 xoffset = detProp.ConvertTicksToX(ticksoffset,allHits[0]->
WireID());
514 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;
520 auto vhit=fmthm.at(i);
521 auto vmeta=fmthm.data(i);
522 for (
size_t ii = 0; ii<vhit.size(); ++ii){
523 bool fBadhit =
false;
529 if (vmeta[ii]->
Index()>=tracklist[i]->NumberTrajectoryPoints()){
530 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!!";
532 if (!tracklist[i]->HasValidPoint(vmeta[ii]->
Index())){
538 auto loc = tracklist[i]->LocationAtPoint(vmeta[ii]->
Index());
544 if (fBadhit)
continue;
545 if (zpos<-100)
continue;
546 planenum=vhit[ii]->WireID().Plane;
548 peakT_0.push_back(vhit[ii]->PeakTime());
551 int_0.push_back(vhit[ii]->Integral());
552 amp_0.push_back(vhit[ii]->PeakAmplitude());
553 rms_0.push_back(vhit[ii]->RMS());
554 hitx_0.push_back(xpos);
555 hity_0.push_back(ypos);
556 hitz_0.push_back(zpos);
559 peakT_1.push_back(vhit[ii]->PeakTime());
562 int_1.push_back(vhit[ii]->Integral());
563 amp_1.push_back(vhit[ii]->PeakAmplitude());
564 rms_1.push_back(vhit[ii]->RMS());
565 hitx_1.push_back(xpos);
566 hity_1.push_back(ypos);
567 hitz_1.push_back(zpos);
571 bool removehit=
false;
572 for (
size_t i2=0; i2<hitlist.size(); ++i2) {
573 if (vhit[ii].
key() == hitlist[i2].key())
continue;
574 if (vhit[ii]->Channel()!=hitlist[i2]->Channel())
continue;
575 if ((vhit[ii]->PeakTime()+30<hitlist[i2]->PeakTime()-30) || (vhit[ii]->PeakTime()-30>hitlist[i2]->PeakTime()+30))
continue;
579 if (removehit)
continue;
582 peakT_2.push_back(vhit[ii]->PeakTime()-this_t0crt2/500.0);
585 int_2.push_back(vhit[ii]->Integral());
586 amp_2.push_back(vhit[ii]->PeakAmplitude());
587 rms_2.push_back(vhit[ii]->RMS());
588 dT_buffer.push_back(vhit[ii]->EndTick()-vhit[ii]->StartTick());
589 hitx_2.push_back(xpos);
590 hity_2.push_back(ypos);
591 hitz_2.push_back(zpos);
592 gof.push_back(vhit[ii]->GoodnessOfFit());
593 multi.push_back(vhit[ii]->Multiplicity());
596 unsigned int wireno=vhit[ii]->WireID().Wire;
598 hitz_wire2.push_back(xyzStart[2]);
602 unsigned int t0=vhit[ii]->PeakTime()-3*(vhit[ii]->RMS());
604 unsigned int t1=vhit[ii]->PeakTime()+3*(vhit[ii]->RMS());
605 if(t1>6000) t1=6000-1;
608 double hit_rms=-9999;
609 std::cout<<
"wirelist size, rawlist size "<<wirelist.size()<<
" "<<rawlist.size()<<
std::endl;
611 for (
unsigned int ich = 0; ich < (rawlist.empty()?wirelist.size():rawlist.size()); ++ich){
615 const auto &wire = wirelist[ich];
616 if(wirelist[ich]->Channel()!=vhit[ii]->Channel())
continue;
618 const auto &signal = wire->Signal();
621 for (
size_t itck = t0; itck <inputsignal.size(); ++itck){
622 if(itck>t1)
continue;
623 inputsignal[itck] = signal[itck];
624 if (inputsignal[itck]>hit_pk){
625 hit_pk = inputsignal[itck];
629 std::cout<<
"hitpeak time, hitraw peak time "<<vhit[ii]->PeakTime()<<
" "<<hit_t<<
std::endl;
631 for(
size_t it1=hit_t-30;it1<hit_t+30;it1++){
632 if(it1<1||it1>5999||it1>inputsignal.size())
continue;
633 hit_signal[ntrks][nhits-1][hitindx]=signal[it1];
640 for (
size_t itck = t0; itck < inputsignal.size(); ++itck){
641 if(itck>t1)
continue;
642 if (inputsignal[itck]>=0.5*hit_pk){
645 if (inputsignal[itck]>=0.1*hit_pk){
646 hit_ch += inputsignal[itck];
647 mean_t += itck*(inputsignal[itck]);
648 mean_t2 += itck*itck*(inputsignal[itck]);
653 hit_rms = sqrt(mean_t2-mean_t*mean_t);
655 else if (!rawlist.empty()){
656 const auto & digitVec = rawlist[ich];
657 if(rawlist[ich]->Channel()!=vhit[ii]->Channel())
continue;
659 raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->GetPedestal(), digitVec->Compression());
662 std::vector<double> signalbuffer;
663 signalbuffer.clear();
665 for(
size_t it=0;it<rawadc.size();it++){
666 signalbuffer.push_back(rawadc[it]);
668 double med_signal=TMath::Median(signalbuffer.size(),&signalbuffer[0]);
669 std::cout<<
"median , pedestal , rawadc size "<<med_signal<<
" "<<digitVec->GetPedestal()<<
" "<<rawadc.size()<<
std::endl;
672 for (
size_t itck = t0; itck <rawadc.size(); ++itck){
673 if(itck>t1)
continue;
675 inputsignal[itck] = rawadc[itck] - med_signal;
676 if (inputsignal[itck]>hit_pk){
677 hit_pk = inputsignal[itck];
681 std::cout<<
"hitpeak time, hitraw peak time "<<vhit[ii]->PeakTime()<<
" "<<hit_t<<
std::endl;
683 for(
size_t it1=hit_t-30;it1<hit_t+30;it1++){
684 if(it1<1|| it1>5999 || it1>rawadc.size())
continue;
686 hit_signal[ntrks][nhits-1][hitindex]=rawadc[it1]-med_signal;
694 for (
size_t itck = t0; itck < inputsignal.size(); ++itck){
695 if(itck>t1)
continue;
697 inputsignal[itck] = rawadc[itck] - med_signal;
698 if (inputsignal[itck]>=0.5*hit_pk){
701 if (inputsignal[itck]>=0.1*hit_pk){
702 hit_ch += inputsignal[itck];
703 mean_t += itck*(inputsignal[itck]);
704 mean_t2 += itck*itck*(inputsignal[itck]);
709 hit_rms = sqrt(mean_t2-mean_t*mean_t);
712 rms_raw2.push_back(hit_rms);
713 hit_peakTrawb.push_back(hit_t);
717 for(
size_t sc=0;sc<fSimChannels.size();++sc){
718 if(fSimChannels[sc]->Channel()==vhit[ii]->Channel()) chan=fSimChannels[sc];
721 const auto &tdcidemap = chan->
TDCIDEMap();
725 double hit_rmsvalue=0;
726 for(
auto mapitr = tdcidemap.begin(); mapitr != tdcidemap.end(); mapitr++){
727 if(!((*mapitr).first>t0 && (*mapitr).first<=t1))
continue;
730 const std::vector<sim::IDE> idevec = (*mapitr).second;
731 for(
size_t iv = 0; iv < idevec.size(); ++iv){
732 hit_nelec+=idevec[iv].numElectrons;
736 int j=(*mapitr).first;
737 mean_t+=double(j)*double(hit_nelec);
738 double jterm=double(j)*double(j)*double(hit_nelec);
739 mean_t2=mean_t2+jterm;
743 hit_rmsvalue=sqrt(mean_t2-mean_t*mean_t);
745 truermsb=hit_rmsvalue;
748 truerms.push_back(truermsb);
755 if(peakT_2.size()<10)
continue;
756 max_value=*std::max_element(peakT_2.begin(),peakT_2.end());
757 min_value=*std::min_element(peakT_2.begin(),peakT_2.end());
759 std::cout<<max_value<<
" "<<min_value<<
std::endl;
807 startcosxyz.push_back(dir_start.X());
808 startcosxyz.push_back(dir_start.Y());
809 startcosxyz.push_back(dir_start.Z());
810 endcosxyz.push_back(dir_end.X());
811 endcosxyz.push_back(dir_end.Y());
812 endcosxyz.push_back(dir_end.Z());
831 t0crt2.push_back(this_t0crt2/500.0);
844 for(
int k=0;
k < 6;
k++)
907 for(
int i=0; i<10; i++){
908 for(
int j=0; j<5000; j++){
909 for(
int k=0;
k<60;
k++){
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 > > trkhitz_wire2
std::vector< float > trkstarty
constexpr std::uint32_t timeLow() const
std::vector< std::vector< float > > hit_rms1
std::vector< std::vector< int > > hit_tpc0
Energy deposited on a readout channel by simulated tracks.
std::vector< std::vector< int > > hit_wire1
std::vector< float > crtreco_y0
hitrmsrawdigits(fhicl::ParameterSet const &pset)
std::vector< std::vector< float > > trkdq_int1
Handle< PROD > getHandle(SelectorBase const &) const
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< std::vector< float > > hit_peakT0
std::vector< std::vector< float > > trkhitx2
std::vector< std::vector< float > > trkhitz2
constexpr std::uint32_t timeHigh() const
std::vector< std::vector< float > > hit_deltaT
unsigned int fWaveformSize
std::vector< float > trackthetaxz
float hit_signal[10][5000][60]
Vector_t VertexDirection() const
std::vector< std::vector< float > > hit_peakT1
std::vector< std::vector< int > > hit_wire2
void analyze(const art::Event &evt)
std::vector< std::vector< float > > hit_peakTraw2
std::vector< float > xprojectedlen
EDAnalyzer(fhicl::ParameterSet const &pset)
std::vector< float > crtreco_y1
std::vector< float > trkendx
std::vector< std::vector< float > > trkhity1
std::vector< float > crtreco_x1
object containing MC flux information
art framework interface to geometry description
std::vector< float > trkendy
std::vector< float > crtreco_x0
std::vector< std::vector< float > > hit_rms_true
std::vector< std::vector< float > > trkdq_amp0
void beginRun(const art::Run &run)
std::vector< std::vector< int > > hit_wire0
std::vector< std::vector< float > > trkhitx1
std::vector< float > trkstartz
double Length(size_t p=0) const
Access to various track properties.
#define DEFINE_ART_MODULE(klass)
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
std::vector< float > trklen
std::vector< float > trkendz
std::vector< std::vector< int > > hit_tpc2
std::vector< std::vector< float > > trkdq_int2
art::InputTag fRawProducerLabel
Point_t const & Vertex() const
std::vector< int > tot_trks
ProtoDUNEDataUtils fDataUtils
SubRunNumber_t subRun() const
std::vector< float > trkstartx_crt2
static int max(int a, int b)
Description of geometry of one entire detector.
std::vector< float > crtreco_z1
Definition of data types for geometry description.
std::vector< std::vector< float > > trkendcosxyz
std::vector< std::vector< float > > multiplicity2
std::vector< std::vector< float > > trkstartcosxyz
std::vector< std::vector< float > > hit_peakT2
std::vector< std::vector< float > > trkdq_amp2
std::vector< std::vector< float > > trkdq_amp1
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
std::string fCalorimetryModuleLabel
std::vector< std::vector< float > > hit_rms0
std::vector< double > crt2tickoffset
Declaration of signal hit object.
std::vector< double > T0_values
art::InputTag fWireProducerLabel
int GetNActiveFembsForAPA(art::Event const &evt, int apa) const
Get number of active fembs in an APA.
std::string fHitsModuleLabel
Vector_t EndDirection() const
std::vector< std::vector< float > > hit_rmsraw2
geo::GeometryCore const * fGeometry
std::vector< std::vector< float > > trkhitx0
Provides recob::Track data product.
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
Declaration of basic channel signal object.
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
virtual ~hitrmsrawdigits()
std::size_t getView(std::string const &moduleLabel, std::string const &productInstanceName, std::string const &processName, std::vector< ELEMENT const * > &result) const
std::vector< std::vector< float > > trkhity2
std::vector< std::vector< float > > hit_rms2
auto const & get(AssnsNode< L, R, D > const &r)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
std::vector< std::vector< float > > trkhitz1
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
recob::tracking::Plane Plane
std::vector< float > trkstartx
std::vector< std::vector< float > > trkhitz0
std::vector< double > t0crt2
std::string fTrackModuleLabel
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
std::vector< float > trkendx_crt2
cet::coded_exception< error, detail::translate > exception
std::vector< std::vector< int > > hit_tpc1
std::vector< float > trackthetayz
QTextStream & endl(QTextStream &s)
Event finding and building.
std::vector< float > crtreco_z0
std::vector< std::vector< float > > trkdq_int0
std::vector< std::vector< float > > goodnessoffit2
std::vector< std::vector< float > > trkhity0