17 #include "TProfile2D.h" 34 #include "Utilities/AssociationUtil.h" 36 #include "Utilities/func/MathUtil.h" 60 #include "TMVA/Tools.h" 61 #include "TMVA/Reader.h" 62 #include "TMVA/MethodCuts.h" 68 namespace rb{
class Cluster;
class Track;}
69 namespace sim{
class Particle;}
70 namespace simb{
class MCFlux;
class MCTruth;
class MCNeutrino;}
217 float prong1cellX[500];
218 float prong1cellY[500];
219 float prong1cellZ[500];
220 float prong1cellE[500];
221 int prong1cellADC[500];
222 int prong1cellPlane[500];
223 int prong1cellCell[500];
224 float prong1cellPECorr[500];
225 float prong1cellW[500];
291 _otree = tfs->
make<TTree>(
"Xeff",
"Xeff particle");
401 _btree=
new TTree(
"beamInfo",
"beam info");
543 for(
int i=0;
i<500; ++
i ){
568 std::vector<int> sq_plane;
570 std::vector<int> sq_xview;
572 std::vector<float> sq_tns;
574 std::vector<int> sq_cell;
577 int sq_ncells=chits->size();
578 for(
int i = 0;
i < sq_ncells; ++
i){
580 double hittime=hit->
TNS()/1000.;
581 sq_plane.push_back(hit->
Plane());
582 sq_tns.push_back(hittime);
583 sq_cell.push_back(hit->
Cell());
584 if( (hit->
View())==
geo::kX ) sq_xview.push_back(1);
585 else sq_xview.push_back(0);
588 int Nmissingdcm[14]={0};
589 for(
int ic=0; ic<sq_ncells; ++ic ){
590 if( sq_plane[ic]<64 ){
591 if( sq_xview[ic]==1 ){
592 if( sq_cell[ic]>63 ) Nmissingdcm[0]++;
593 else Nmissingdcm[1]++;
596 if( sq_cell[ic]>31 ) Nmissingdcm[2]++;
597 else Nmissingdcm[3]++;
600 else if( sq_plane[ic]<128 ){
601 if( sq_xview[ic]==1 ){
602 if( sq_cell[ic]>63 ) Nmissingdcm[4]++;
603 else Nmissingdcm[5]++;
606 if( sq_cell[ic]>31 ) Nmissingdcm[6]++;
607 else Nmissingdcm[7]++;
610 else if( sq_plane[ic]<192 ){
611 if( sq_xview[ic]==1 ){
612 if( sq_cell[ic]>63 ) Nmissingdcm[8]++;
613 else Nmissingdcm[9]++;
616 if( sq_cell[ic]>31 ) Nmissingdcm[10]++;
617 else Nmissingdcm[11]++;
621 if( sq_xview[ic]==1 ) Nmissingdcm[12]++;
622 else Nmissingdcm[13]++;
627 for(
int id=0;
id<14; ++
id ){
628 if( Nmissingdcm[
id]==0 ) ndeaddcms++;
649 for(
unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
662 for(
unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
670 for(
size_t h = 0;
h < chits->size(); ++
h)
687 std::vector< art::Ptr<simb::MCFlux> > fluxs = fmf.at(0);
689 if(fluxs.size() == 1){
699 if( !(nu.
CCNC()==0) )
continue;
721 if(
fabs(nu.
Nu().
Vx())>200. )
continue;
722 if(
fabs(nu.
Nu().
Vy())>200. )
continue;
723 if( nu.
Nu().
Vz()<0. )
continue;
735 if( mclist->size()>0 && fllist->size()>0 ){
736 for(
unsigned int i_intx=0; i_intx<mclist->size(); ++i_intx){
739 TLorentzVector Tslnu(nu.
Nu().
Px(),nu.
Nu().
Py(),
741 TLorentzVector Tmcnu(mcnu.
Nu().
Px(),mcnu.
Nu().
Py(),
762 if (!(fmVertex.isValid())){
767 std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
768 if( vert.size() != 1 ){
776 std::vector<art::Ptr<rb::Track>> Ktracks = fmKTrack.at(sliceIdx);
777 double longestKTrack = 0;
778 int longestKTrack_ncells=0;
779 double longestKTrack_energy=0.;
780 int longestKTrack_nplanes=0;
781 float longestKTrack_phi=0.;
782 float longestKTrack_theta=0.;
783 float longestKTrack_startx=0.;
784 float longestKTrack_starty=0.;
785 float longestKTrack_startz=0.;
786 float longestKTrack_stopx=0.;
787 float longestKTrack_stopy=0.;
788 float longestKTrack_stopz=0.;
789 int longestKTrack_is3D=0;
791 if( fmKTrack.isValid() ){
792 for(
unsigned int n = 0;
n < Ktracks.size(); ++
n){
796 const double L = Ktracks[
n]->TotalLength();
797 if(L > longestKTrack){
799 longestKTrack_ncells=Ktracks[
n]->NCell();
800 longestKTrack_energy=Ktracks[
n]->TotalGeV();
801 longestKTrack_nplanes=Ktracks[
n]->ExtentPlane();
802 TVector3
dir = Ktracks[
n]->Dir();
803 longestKTrack_phi=dir.Phi();
804 longestKTrack_theta=dir.Theta();
805 TVector3
start = Ktracks[
n]->Start();
806 TVector3 stop = Ktracks[
n]->Stop();
807 longestKTrack_startx = start.X();
808 longestKTrack_starty = start.Y();
809 longestKTrack_startz = start.Z();
810 longestKTrack_stopx = stop.X();
811 longestKTrack_stopy = stop.Y();
812 longestKTrack_stopz = stop.Z();
813 if( Ktracks[
n]->Is3D() ) longestKTrack_is3D=1;
819 double nCells = slice.
NCell();
823 double recoEnergy=0.;
826 for(
unsigned int icell=0; icell<slice.
NCell(); ++ icell){
830 if( !rhit.IsCalibrated() )
continue;
832 recoEnergy += rhit.
GeV();
834 if( rhit.PECorr()>100. && rhit.PECorr()<245. ) nmip++;
875 std::vector<art::Ptr<rb::Prong>> SelectedProng = fmProng3D.at(0);
877 int nprongs=SelectedProng.size();
878 if( nprongs>20 ) nprongs=20;
880 std::vector<double> ProngEnergy;
881 std::vector<double> ProngEnergyXView;
882 std::vector<double> ProngEnergyYView;
883 std::vector<double> ProngNCells;
884 std::vector<double> ProngNCellsXView;
885 std::vector<double> ProngNCellsYView;
886 std::vector<double> ProngPhi;
887 std::vector<double> ProngTheta;
888 std::vector<double> ProngLength;
889 std::vector<double> ProngStartX;
890 std::vector<double> ProngStartY;
891 std::vector<double> ProngStartZ;
892 std::vector<double> ProngStopX;
893 std::vector<double> ProngStopY;
894 std::vector<double> ProngStopZ;
895 std::vector<int> ProngNplanes;
903 if( SelectedProng.size()>0 ){
905 std::vector<int> prongIndex;
906 std::vector<float> pcellPECorr;
910 TVector3 dir1 = SelectedProng[
ip]->Dir();
911 TVector3 start1 = SelectedProng[
ip]->Start();
912 TVector3 Pstop = start1 + (SelectedProng[
ip]->TotalLength())*dir1;
923 for(
unsigned int icell=0; icell<SelectedProng[
ip]->NCell(); ++ icell){
927 if( !rhit.IsCalibrated() )
continue;
928 Eprong += rhit.
GeV();
932 EprongX += rhit.GeV();
936 EprongY += rhit.GeV();
939 prongIndex.push_back(
ip);
940 pcellPECorr.push_back(rhit.PECorr());
942 ProngTheta.push_back(dir1.Theta());
943 ProngLength.push_back(SelectedProng[
ip]->TotalLength());
945 ProngStartX.push_back(start1.X());
946 ProngStartY.push_back(start1.Y());
947 ProngStartZ.push_back(start1.Z());
948 ProngStopX.push_back(Pstop.X());
949 ProngStopY.push_back(Pstop.Y());
950 ProngStopZ.push_back(Pstop.Z());
952 ProngEnergy.push_back(Eprong);
953 ProngEnergyXView.push_back(EprongX);
954 ProngEnergyYView.push_back(EprongY);
955 ProngPhi.push_back(dir1.Phi());
956 ProngNCells.push_back(npcells);
957 ProngNCellsXView.push_back(nxpcells);
958 ProngNCellsYView.push_back(nypcells);
960 ProngNplanes.push_back(SelectedProng[
ip]->ExtentPlane());
964 double MaxEnergy1=-999.;
965 double MaxEnergy2=-999.;
967 if( MaxEnergy1<ProngEnergy[
i] ){
968 MaxEnergy1=ProngEnergy[
i];
973 p1Phi=ProngPhi[Index1];
993 std::vector<float> p1cellX;
994 std::vector<float> p1cellY;
995 std::vector<float> p1cellZ;
996 std::vector<float> p1cellE;
997 std::vector<int> p1cellPlane;
998 std::vector<int> p1cellCell;
999 std::vector<int> p1cellADC;
1000 std::vector<float> p1cellPE;
1001 std::vector<float> p1cellPECorr;
1002 std::vector<int> p1cellXView;
1003 std::vector<float> p1cellTNS;
1004 std::vector<float> p1cellW;
1006 for(
unsigned int icell=0; icell<SelectedProng[Index1]->NCell(); ++ icell){
1009 double wx=SelectedProng[Index1]->W(chit.
get());
1011 if( !rhit.IsCalibrated() )
continue;
1013 if( rhit.PECorr()>100. && rhit.PECorr()<240. ) ++p1_mipcells;
1016 p1cellW.push_back(wx);
1019 p1cellW.push_back(wy);
1022 p1cellX.push_back(rhit.X());
1023 p1cellY.push_back(rhit.Y());
1024 p1cellZ.push_back(rhit.Z());
1025 p1cellE.push_back(rhit.GeV());
1027 p1cellADC.push_back(chit->
ADC());
1028 p1cellPlane.push_back(chit->
Plane());
1029 p1cellCell.push_back(chit->
Cell());
1030 p1cellPECorr.push_back(rhit.PECorr());
1037 if( p1_ncells>500 ) p1_ncells=500;
1040 for(
int ic=0; ic<p1_ncells; ++ic ){
1052 if( minplane>p1cellPlane[ic] ){
1053 minplane=p1cellPlane[ic];
1056 if( p1_ncells>0 )
p1Fmip=(
float)p1_mipcells/p1_ncells;
1058 int Nplanes_p1=ProngNplanes[Index1];
1060 if( ProngEnergy[Index1]>0. ){
1061 std::vector<double> p1PlaneEnergy;
1062 p1PlaneEnergy.clear();
1063 double totalE_road20=0.;
1064 for(
int ip=0;
ip<Nplanes_p1; ++
ip ){
1073 std::vector<double> Xxview;
1074 std::vector<double> Exview;
1075 std::vector<double> Yyview;
1076 std::vector<double> Eyview;
1082 for(
int ic=0; ic<p1_ncells; ++ic ){
1083 int Dplane=p1cellPlane[ic]-minplane;
1085 eplane += p1cellE[ic];
1087 Ecell += p1cellE[ic];
1088 EXcell += p1cellE[ic]*p1cellX[ic];
1089 EYcell += p1cellE[ic]*p1cellY[ic];
1091 if( p1cellPlane[ic]%2 != 0 ){
1092 Xxview.push_back(p1cellX[ic]);
1093 Exview.push_back(p1cellE[ic]);
1096 Yyview.push_back(p1cellY[ic]);
1097 Eyview.push_back(p1cellE[ic]);
1102 p1PlaneEnergy.push_back(eplane);
1105 cellX = EXcell/Ecell;
1106 cellY = EYcell/Ecell;
1108 int Nxviews=Xxview.size();
1109 int Nyviews=Yyview.size();
1111 for(
int ixv=0; ixv<Nxviews; ++ixv ){
1112 if(
fabs(Xxview[ixv]-cellX) < 2.*2. ) totalE_road20+=Exview[ixv];
1116 for(
int iyv=0; iyv<Nyviews; ++iyv ){
1117 if(
fabs(Yyview[iyv]-cellY) < 2.*2. ) totalE_road20+=Eyview[iyv];
1123 double efrac_6plane=-1.;
1124 double p1_Eplane[20]={0.};
1126 int nplanes=p1PlaneEnergy.size();
1128 for(
int j=0;
j<20; ++
j ){
1129 if(
ip<=
j ) p1_Eplane[
j]+=p1PlaneEnergy[
ip];
1133 double e_frac6=(p1PlaneEnergy[
ip]+p1PlaneEnergy[
ip+1]+p1PlaneEnergy[
ip+2]+p1PlaneEnergy[
ip+3]+p1PlaneEnergy[
ip+4]+p1PlaneEnergy[
ip+5])/ProngEnergy[Index1];
1134 if( efrac_6plane<e_frac6 ) efrac_6plane=e_frac6;
1136 else if( nplanes<=6 ) efrac_6plane=1.;
1140 double efrac_Eplane[20]={0.};
1141 for(
int ip=0;
ip<20; ++
ip ){
1142 efrac_Eplane[
ip] = p1_Eplane[
ip]/ProngEnergy[Index1];
1145 double efrac_plane2=0.;
1146 double efrac_plane3=0.;
1147 double efrac_plane4=0.;
1149 if( p1PlaneEnergy.size()>1 && ProngEnergy[Index1]>0. )
1150 efrac_plane2=p1PlaneEnergy[1]/ProngEnergy[Index1];
1151 if( p1PlaneEnergy.size()>2 && ProngEnergy[Index1]>0. )
1152 efrac_plane3=p1PlaneEnergy[2]/ProngEnergy[Index1];
1153 if( p1PlaneEnergy.size()>3 && ProngEnergy[Index1]>0. )
1154 efrac_plane4=p1PlaneEnergy[3]/ProngEnergy[Index1];
1161 efrac_2sig=totalE_road20/ProngEnergy[Index1];
1168 if(
i==Index1 )
continue;
1169 if( MaxEnergy2<ProngEnergy[
i] ){
1170 MaxEnergy2=ProngEnergy[
i];
1174 double Eden=ProngEnergy[Index1]+ProngEnergy[Index2];
1175 double Enum=ProngEnergy[Index1]-ProngEnergy[Index2];
1176 if( Eden>0.0 ) ratio=Enum/Eden;
1177 dphi=
fabs(ProngPhi[Index1]-ProngPhi[Index2]);
1178 if( dphi>2.*TMath::Pi() ) dphi=dphi-2.*TMath::Pi();
1179 if( dphi>TMath::Pi() ) dphi=2.*TMath::Pi()-dphi;
1182 if( prongIndex.size()>0 ){
1184 double ntotal_p2=0.;
1185 for(
unsigned int ic=0; ic<prongIndex.size(); ++ic ){
1186 if( prongIndex[ic] != Index2 )
continue;
1188 if( pcellPECorr[ic]>100. && pcellPECorr[ic]<240. )
1191 if( nmips_p2>0. )
p2Fmip = nmips_p2/ntotal_p2;
1205 lid = elid->Value();
1220 TH1* hPOT = tfs->
make<TH1F>(
"hTotalPOT",
";; POT", 1, 0, 1);
1222 TH1* hSPILL = tfs->
make<TH1F>(
"hTotalSpill",
";; SPILL", 1, 0, 1);
std::string fProngLabel
label for kalman tracks
double E(const int i=0) const
::xsd::cxx::tree::id< char, ncname > id
SubRunNumber_t subRun() const
back track the reconstruction to the simulation
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
double Theta() const
angle between incoming and outgoing leptons, in radians
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
const simb::MCNeutrino & GetNeutrino() const
double Py(const int i=0) const
float prong1cellPECorr[500]
fvar< T > fabs(const fvar< T > &x)
SubRunNumber_t subRun() const
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned short Plane() const
std::string fKTrackLabel
label for slices
void analyze(const art::Event &evt)
const simb::MCParticle & Nu() const
Vertical planes which measure X.
double Pt() const
transverse momentum of interaction, in GeV/c
double Px(const int i=0) const
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
std::string fCellHitLabel
label for vertex
charged current quasi-elastic
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
Horizontal planes which measure Y.
object containing MC flux information
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
void reconfigure(const fhicl::ParameterSet &p)
std::string fInstLabel2D
instance label for prongs 3D
std::string fMCFluxLabel
instance label for prongs 2D
unsigned short Cell() const
int InteractionType() const
DEFINE_ART_MODULE(ROCKMRE)
const simb::MCParticle & Lepton() const
void push_back(Ptr< U > const &p)
void endSubRun(const art::SubRun &sr)
static constexpr double L
T get(std::string const &key) const
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
const rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
This class describes a particle created in the detector Monte Carlo simulation.
std::string fVertexLabel
label for generator information
EventNumber_t event() const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
double Vx(const int i=0) const
T * make(ARGS...args) const
int16_t ADC(uint32_t i) const
Data resulting from a Hough transform on the cell hit positions.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
double Pz(const int i=0) const
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
double Vz(const int i=0) const
std::string fInstLabel3D
label for prong
bool IsNoise() const
Is the noise flag set?
double totgoodpot
normalized by 10^12 POT
Event generator information.
double Vy(const int i=0) const
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string flidLabel
label for cell hits