17 #include "TProfile2D.h" 34 #include "Utilities/AssociationUtil.h" 36 #include "Utilities/func/MathUtil.h" 61 #include "TMVA/Tools.h" 62 #include "TMVA/Reader.h" 63 #include "TMVA/MethodCuts.h" 69 namespace rb{
class Cluster;
class Track;}
70 namespace sim{
class Particle;}
71 namespace simb{
class MCFlux;
class MCTruth;
class MCNeutrino;}
145 float PlaneEnergy[192];
147 int PlaneNcells[192];
161 int slcellPlane[500];
163 float slcellPECorr[500];
171 float prongEnergy[20];
173 float prongTheta[20];
174 float prongStartX[20];
175 float prongStartY[20];
176 float prongStartZ[20];
177 float prongLength[20];
178 float prongStopX[20];
179 float prongStopY[20];
180 float prongStopZ[20];
181 float prongEnergyXView[20];
182 float prongEnergyYView[20];
184 int prongNCellsXView[20];
185 int prongNCellsYView[20];
186 int prongNplanes[20];
191 float prong1cellX[500];
192 float prong1cellY[500];
193 float prong1cellZ[500];
194 float prong1cellE[500];
195 int prong1cellADC[500];
196 int prong1cellPlane[500];
197 int prong1cellCell[500];
198 float prong1cellPECorr[500];
199 float prong1cellW[500];
203 int prong2DInXView[20];
204 float prongEnergy2D[20];
205 int prongNCells2D[20];
206 float prongPhi2D[20];
207 float prongTheta2D[20];
208 float prongStartX2D[20];
209 float prongStartY2D[20];
210 float prongStartZ2D[20];
211 float prongLength2D[20];
212 float prongStopX2D[20];
213 float prongStopY2D[20];
214 float prongStopZ2D[20];
248 isMC = pset.
get<
int>(
"MCevts");
260 _otree = tfs->
make<TTree>(
"ROCKMRE",
"ROCKMRE particle");
355 _btree=
new TTree(
"beamInfo",
"beam info");
424 for(
int i=0;
i<192; ++
i ){
436 for(
int i=0;
i<20; ++
i ){
473 for(
int i=0;
i<500; ++
i ){
519 std::vector<int> sq_plane;
521 std::vector<int> sq_xview;
523 std::vector<float> sq_tns;
525 std::vector<int> sq_cell;
528 int sq_ncells=chits0->size();
529 for(
int i = 0;
i < sq_ncells; ++
i){
531 double hittime=hit->
TNS()/1000.;
532 sq_plane.push_back(hit->
Plane());
533 sq_tns.push_back(hittime);
534 sq_cell.push_back(hit->
Cell());
535 if( (hit->
View())==
geo::kX ) sq_xview.push_back(1);
536 else sq_xview.push_back(0);
540 int Ndbdcm_intime[14]={0};
541 int Ndbdcm_outtime[14]={0};
543 for(
int ic=0; ic<sq_ncells; ++ic ){
544 if( sq_plane[ic]<64 ){
545 if( sq_xview[ic]==1 ){
546 if( sq_cell[ic]>63 ){
547 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[0]++;
548 else Ndbdcm_outtime[0]++;
551 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[1]++;
552 else Ndbdcm_outtime[1]++;
556 if( sq_cell[ic]>31 ){
557 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[2]++;
558 else Ndbdcm_outtime[2]++;
561 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[3]++;
562 else Ndbdcm_outtime[3]++;
566 else if( sq_plane[ic]<128 ){
567 if( sq_xview[ic]==1 ){
568 if( sq_cell[ic]>63 ){
569 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[4]++;
570 else Ndbdcm_outtime[4]++;
573 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[5]++;
574 else Ndbdcm_outtime[5]++;
578 if( sq_cell[ic]>31 ){
579 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[6]++;
580 else Ndbdcm_outtime[6]++;
583 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[7]++;
584 else Ndbdcm_outtime[7]++;
588 else if( sq_plane[ic]<192 ){
589 if( sq_xview[ic]==1 ){
590 if( sq_cell[ic]>63 ){
591 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[8]++;
592 else Ndbdcm_outtime[8]++;
595 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[9]++;
596 else Ndbdcm_outtime[9]++;
600 if( sq_cell[ic]>31 ){
601 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[10]++;
602 else Ndbdcm_outtime[10]++;
605 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[11]++;
606 else Ndbdcm_outtime[11]++;
611 if( sq_xview[ic]==1 ){
612 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[12]++;
613 else Ndbdcm_outtime[12]++;
616 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[13]++;
617 else Ndbdcm_outtime[13]++;
623 std::vector<std::string> dcmlocation;
625 for(
int id=0;
id<14; ++
id ){
626 if( Ndbdcm_outtime[
id]==0 && Ndbdcm_intime[
id]==0 ){
639 for(
unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
654 for(
unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
659 if (!(fmVertex.isValid()))
continue;
660 std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
662 if (vert.size() != 1)
continue;
666 double nCells = slice.
NCell();
670 double recoEnergy=0.;
673 for(
unsigned int icell=0; icell<slice.
NCell(); ++ icell){
677 if( !rhit.IsCalibrated() )
continue;
678 recoEnergy += rhit.
GeV();
680 if( rhit.PECorr()>100. && rhit.PECorr()<245. ) nmip++;
698 if( nCells<10 || nCells>500 )
continue;
761 std::vector<art::Ptr<rb::Prong>> SelectedProng = fmProng3D.at(0);
764 if( SelectedProng.size()<1 )
continue;
765 int nprongs=SelectedProng.size();
766 if( nprongs>20 ) nprongs=20;
768 std::vector<double> ProngEnergy;
769 std::vector<double> ProngEnergyXView;
770 std::vector<double> ProngEnergyYView;
771 std::vector<double> ProngPhi;
776 TVector3 dir1 = SelectedProng[
ip]->Dir();
781 TVector3 start1 = SelectedProng[
ip]->Start();
788 TVector3 Pstop = start1 + (SelectedProng[
ip]->TotalLength())*dir1;
804 for(
unsigned int icell=0; icell<SelectedProng[
ip]->NCell(); ++ icell){
808 if( !rhit.IsCalibrated() )
continue;
809 Eprong += rhit.
GeV();
813 EprongX += rhit.GeV();
817 EprongY += rhit.GeV();
822 ProngEnergy.push_back(Eprong);
823 ProngEnergyXView.push_back(EprongX);
824 ProngEnergyYView.push_back(EprongY);
825 ProngPhi.push_back(dir1.Phi());
835 double MaxEnergy1=-999.;
838 if( MaxEnergy1<ProngEnergy[
i] ){
839 MaxEnergy1=ProngEnergy[
i];
849 std::vector<float> p1cellX;
850 std::vector<float> p1cellY;
851 std::vector<float> p1cellZ;
852 std::vector<float> p1cellE;
853 std::vector<int> p1cellPlane;
854 std::vector<int> p1cellCell;
855 std::vector<int> p1cellADC;
856 std::vector<float> p1cellPE;
857 std::vector<float> p1cellPECorr;
858 std::vector<int> p1cellXView;
859 std::vector<float> p1cellTNS;
860 std::vector<float> p1cellW;
862 for(
unsigned int icell=0; icell<SelectedProng[Index1]->NCell(); ++ icell){
865 double wx=SelectedProng[Index1]->W(chit.
get());
867 if( !rhit.IsCalibrated() )
continue;
869 if( rhit.PECorr()>100. && rhit.PECorr()<240. ) ++p1_mipcells;
872 p1cellW.push_back(wx);
875 p1cellW.push_back(wy);
878 p1cellX.push_back(rhit.X());
879 p1cellY.push_back(rhit.Y());
880 p1cellZ.push_back(rhit.Z());
881 p1cellE.push_back(rhit.GeV());
883 p1cellADC.push_back(chit->
ADC());
884 p1cellPlane.push_back(chit->
Plane());
885 p1cellCell.push_back(chit->
Cell());
886 p1cellPECorr.push_back(rhit.PECorr());
891 for(
int ic=0; ic<p1_ncells; ++ic ){
907 std::vector<float> slicecellX;
908 std::vector<float> slicecellY;
909 std::vector<float> slicecellZ;
910 std::vector<float> slicecellE;
911 std::vector<int> slicecellPlane;
912 std::vector<int> slicecellCell;
913 std::vector<int> slicecellADC;
914 std::vector<float> slicecellPE;
915 std::vector<float> slicecellPECorr;
916 std::vector<int> slicecellXView;
917 std::vector<float> slicecellTNS;
918 std::vector<float> slicecellW;
919 std::vector<int> slicecellNoise;
921 for(
unsigned int icell=0; icell<slice.
NCell(); ++ icell){
926 if( !rhit.IsCalibrated() )
continue;
930 slicecellW.push_back(slice.
W(chit.
get()));
932 slicecellX.push_back(rhit.X());
933 slicecellY.push_back(rhit.Y());
934 slicecellZ.push_back(rhit.Z());
935 slicecellE.push_back(rhit.GeV());
937 slicecellADC.push_back(chit->
ADC());
938 slicecellPlane.push_back(chit->
Plane());
939 slicecellCell.push_back(chit->
Cell());
940 slicecellPECorr.push_back(rhit.PECorr());
946 for(
int ic=0; ic<ncells_slice; ++ic ){
963 std::vector<art::Ptr<rb::Prong>> SelectedProng2D = fmProng2D.at(0);
965 std::vector<int> Prong2DInXView;
966 std::vector<double> ProngEnergy2D;
967 int nprongs2D=SelectedProng2D.size();
968 if( nprongs2D>20 ) nprongs2D=20;
971 for(
int ip=0;
ip<nprongs2D; ++
ip ){
973 int prongview=SelectedProng2D[
ip]->View();
978 TVector3 dir1 = SelectedProng2D[
ip]->Dir();
982 TVector3 start1 = SelectedProng2D[
ip]->Start();
987 TVector3 Pstop = start1 + (SelectedProng2D[
ip]->TotalLength())*dir1;
992 Prong2DInXView.push_back(prongview);
995 for(
unsigned int icell=0; icell<SelectedProng2D[
ip]->NCell(); ++ icell){
997 const TVector3
mean = SelectedProng2D[
ip]->MeanXYZ();
1000 if( !rhit.IsCalibrated() )
continue;
1001 Eprong += rhit.
GeV();
1005 ProngEnergy2D.push_back(Eprong);
1011 if( max3D<ProngEnergy[
i] ){
1012 max3D=ProngEnergy[
i];
1018 for(
int i=0;
i<nprongs2D; ++
i ){
1019 if( max2D<ProngEnergy2D[
i] ){
1020 max2D=ProngEnergy2D[
i];
1024 double ebalance2D=1.;
1025 if( nprongs>0 && nprongs2D>0 ){
1028 if( Prong2DInXView[index2D]==0 ){
1029 Eden=ProngEnergy2D[index2D]+ProngEnergyXView[index3D];
1030 Enum=ProngEnergy2D[index2D]-ProngEnergyXView[index3D];
1033 Eden=ProngEnergy2D[index2D]+ProngEnergyYView[index3D];
1034 Enum=ProngEnergy2D[index2D]-ProngEnergyYView[index3D];
1036 if( Eden>0. ) ebalance2D=
fabs(Enum)/Eden;
1053 TH1* hPOT = tfs->
make<TH1F>(
"hTotalPOT",
";; POT", 1, 0, 1);
1055 TH1* hSPILL = tfs->
make<TH1F>(
"hTotalSpill",
";; SPILL", 1, 0, 1);
::xsd::cxx::tree::id< char, ncname > id
void analyze(const art::Event &evt)
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.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
fvar< T > fabs(const fvar< T > &x)
SubRunNumber_t subRun() const
float prong1cellPECorr[500]
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned short Plane() const
float prongEnergyYView[20]
Vertical planes which measure X.
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
std::string fInstLabel3D
label for prong
std::string fProngLabel
label for slices
const rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
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...
unsigned short Cell() const
float prongEnergyXView[20]
DEFINE_ART_MODULE(ROCKMRE)
std::string fCellHitLabel
label for vertex
T get(std::string const &key) const
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
void reconfigure(const fhicl::ParameterSet &p)
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
This class describes a particle created in the detector Monte Carlo simulation.
std::string fVertexLabel
label for generator information
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
EventNumber_t event() const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
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
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
std::string fDataSpillLabel
label for cell hits
bool IsNoise() const
Is the noise flag set?
double totgoodpot
normalized by 10^12 POT
std::string fInstLabel2D
instance label for prongs 3D
void endSubRun(const art::SubRun &sr)
double spillpot
POT for spill normalized by 10^12.
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string fMCFluxLabel
instance label for prongs 2D