21 #include "TProfile2D.h" 40 #include "Utilities/AssociationUtil.h" 42 #include "Utilities/func/MathUtil.h" 68 #include "TMVA/Tools.h" 69 #include "TMVA/Reader.h" 70 #include "TMVA/MethodCuts.h" 76 namespace rb{
class Cluster;
class Track;}
77 namespace sim{
class Particle;}
78 namespace simb{
class MCFlux;
class MCTruth;
class MCNeutrino;}
223 float plus1sigma[68];
224 float plus2sigma[68];
225 float minus1sigma[68];
226 float minus2sigma[68];
254 float PlaneEnergy[192];
256 int PlaneNcells[192];
281 float trackcellX[500];
282 float trackcellY[500];
283 float trackcellZ[500];
284 float trackcellE[500];
285 int trackcellADC[500];
286 float trackcellTNS[500];
287 int trackcellPlane[500];
288 int trackcellCell[500];
289 int trackcellView[500];
290 float trackcellPE[500];
291 float trackcellPECorr[500];
292 float trackcellW[500];
305 int slcellPlane[500];
307 float slcellPECorr[500];
309 float slcellTNS[500];
316 float prongEnergy[20];
318 float prongTheta[20];
319 float prongStartX[20];
320 float prongStartY[20];
321 float prongStartZ[20];
322 float prongLength[20];
323 float prongStopX[20];
324 float prongStopY[20];
325 float prongStopZ[20];
326 float prongEnergyXView[20];
327 float prongEnergyYView[20];
329 int prongNCellsXView[20];
330 int prongNCellsYView[20];
331 int prongNplanes[20];
334 int allprongIndex[500];
335 float allprongcellX[500];
336 float allprongcellY[500];
337 float allprongcellZ[500];
338 float allprongcellE[500];
339 int allprongcellADC[500];
340 int allprongcellPlane[500];
341 int allprongcellCell[500];
342 float allprongcellPECorr[500];
343 float allprongcellW[500];
371 float prong1cellX[500];
372 float prong1cellY[500];
373 float prong1cellZ[500];
374 float prong1cellE[500];
375 int prong1cellADC[500];
376 int prong1cellPlane[500];
377 int prong1cellCell[500];
378 float prong1cellPECorr[500];
379 float prong1cellW[500];
400 int prong2DInXView[20];
401 float prongEnergy2D[20];
402 int prongNCells2D[20];
403 float prongPhi2D[20];
404 float prongTheta2D[20];
405 float prongStartX2D[20];
406 float prongStartY2D[20];
407 float prongStartZ2D[20];
408 float prongLength2D[20];
409 float prongStopX2D[20];
410 float prongStopY2D[20];
411 float prongStopZ2D[20];
415 int allprong2DIndex[500];
416 float allprong2DcellX[500];
417 float allprong2DcellY[500];
418 float allprong2DcellZ[500];
419 float allprong2DcellE[500];
420 int allprong2DcellADC[500];
421 int allprong2DcellPlane[500];
422 int allprong2DcellCell[500];
423 float allprong2DcellPECorr[500];
433 int shwNCellsXview[20];
434 int shwNCellsYview[20];
453 int shwcellPlane[500];
454 int shwcellCell[500];
455 float shwcellPECorr[500];
507 isMC = pset.
get<
int>(
"MCevts");
530 _otree = tfs->
make<TTree>(
"NCAna",
"NCAna particle");
793 _otree->Branch(
"shw3D",&
shw3D,
"shw3D[Nshowers]/I");
808 _btree=
new TTree(
"beamInfo",
"beam info");
898 for(
int i=0;
i<14; ++
i )
926 for(
int i=0;
i<14; ++
i )
943 for(
int i=0;
i<20; ++
i ){
993 for(
int i=0;
i<68; ++
i ){
1012 for(
int i=0;
i<20; ++
i ){
1030 for(
int i=0;
i<192; ++
i ){
1058 for(
int i=0;
i<20; ++
i ){
1102 for(
int i=0;
i<500; ++
i ){
1199 std::string dcmName[14]={
"DB1-DCM1",
"DB1-DCM2",
"DB1-DCM3",
"DB1-DCM4",
1200 "DB2-DCM1",
"DB2-DCM2",
"DB2-DCM3",
"DB2-DCM4",
1201 "DB3-DCM1",
"DB3-DCM2",
"DB3-DCM3",
"DB3-DCM4",
1202 "MC-DCM1",
"MC-DCM2"};
1236 std::vector<int> sq_plane;
1238 std::vector<int> sq_xview;
1240 std::vector<float> sq_tns;
1242 std::vector<int> sq_cell;
1245 int sq_ncells=chits->size();
1246 for(
int i = 0;
i < sq_ncells; ++
i){
1248 double hittime=hit->
TNS()/1000.;
1249 sq_plane.push_back(hit->
Plane());
1250 sq_tns.push_back(hittime);
1251 sq_cell.push_back(hit->
Cell());
1252 if( (hit->
View())==
geo::kX ) sq_xview.push_back(1);
1253 else sq_xview.push_back(0);
1256 int ncells_intime=0;
1257 int ncells_outtime=0;
1259 int Ndbdcm_intime[14]={0};
1260 int Ndbdcm_outtime[14]={0};
1262 int Ndb1apd[192]={0};
1263 int Ndb2apd[192]={0};
1264 int Ndb3apd[192]={0};
1267 for(
int ic=0; ic<sq_ncells; ++ic ){
1268 if( sq_tns[ic]>218. && sq_tns[ic]<228. ){
1274 if( sq_plane[ic]<64 ){
1275 for(
int ip=0;
ip<64; ++
ip ){
1276 if(
ip == sq_plane[ic] ){
1277 if( sq_cell[ic]<32 ) Ndb1apd[0+
ip*3]++;
1278 else if( sq_cell[ic]<64 ) Ndb1apd[1+
ip*3]++;
1279 else Ndb1apd[2+
ip*3]++;
1283 else if( sq_plane[ic]<128 ){
1284 for(
int ip=0;
ip<64; ++
ip ){
1285 if( (
ip+64) == sq_plane[ic] ){
1286 if( sq_cell[ic]<32 ) Ndb2apd[0+
ip*3]++;
1287 else if( sq_cell[ic]<64 ) Ndb2apd[1+
ip*3]++;
1288 else Ndb2apd[2+
ip*3]++;
1292 else if( sq_plane[ic]<192 ){
1293 for(
int ip=0;
ip<64; ++
ip ){
1294 if( (
ip+128) == sq_plane[ic] ){
1295 if( sq_cell[ic]<32 ) Ndb3apd[0+
ip*3]++;
1296 else if( sq_cell[ic]<64 ) Ndb3apd[1+
ip*3]++;
1297 else Ndb3apd[2+
ip*3]++;
1302 for(
int ip=0;
ip<22; ++
ip ){
1303 if( (
ip+192) == sq_plane[ic] ){
1304 if( (
ip+192)%2 ==0 ){
1305 if( sq_cell[ic]<32 ) Nmcapd[0+5*
ip/2]++;
1306 else Nmcapd[1+5*
ip/2]++;
1309 if( sq_cell[ic]<32 ) Nmcapd[2+5*(
ip-1)/2]++;
1310 else if( sq_cell[ic]<64 ) Nmcapd[3+5*(
ip-1)/2]++;
1311 else Nmcapd[4+5*(
ip-1)/2]++;
1317 if( sq_plane[ic]<64 ){
1318 if( sq_xview[ic]==1 ){
1319 if( sq_cell[ic]>63 ){
1320 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[0]++;
1321 else Ndbdcm_outtime[0]++;
1324 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[1]++;
1325 else Ndbdcm_outtime[1]++;
1329 if( sq_cell[ic]>31 ){
1330 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[2]++;
1331 else Ndbdcm_outtime[2]++;
1334 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[3]++;
1335 else Ndbdcm_outtime[3]++;
1339 else if( sq_plane[ic]<128 ){
1340 if( sq_xview[ic]==1 ){
1341 if( sq_cell[ic]>63 ){
1342 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[4]++;
1343 else Ndbdcm_outtime[4]++;
1346 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[5]++;
1347 else Ndbdcm_outtime[5]++;
1351 if( sq_cell[ic]>31 ){
1352 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[6]++;
1353 else Ndbdcm_outtime[6]++;
1356 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[7]++;
1357 else Ndbdcm_outtime[7]++;
1361 else if( sq_plane[ic]<192 ){
1362 if( sq_xview[ic]==1 ){
1363 if( sq_cell[ic]>63 ){
1364 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[8]++;
1365 else Ndbdcm_outtime[8]++;
1368 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[9]++;
1369 else Ndbdcm_outtime[9]++;
1373 if( sq_cell[ic]>31 ){
1374 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[10]++;
1375 else Ndbdcm_outtime[10]++;
1378 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[11]++;
1379 else Ndbdcm_outtime[11]++;
1384 if( sq_xview[ic]==1 ){
1385 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[12]++;
1386 else Ndbdcm_outtime[12]++;
1389 if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[13]++;
1390 else Ndbdcm_outtime[13]++;
1396 std::vector<std::string> dcmlocation;
1397 dcmlocation.clear();
1398 for(
int id=0;
id<14; ++
id ){
1399 if( Ndbdcm_outtime[
id]==0 && Ndbdcm_intime[
id]==0 ){
1401 dcmlocation.push_back(dcmName[
id]);
1405 int Napd_24cells[4]={0};
1406 int Napd_28cells[4]={0};
1407 for(
int ia=0; ia<192; ++ia ){
1409 if( Ndb1apd[ia]>=24 ) Napd_24cells[0]+=1.;
1410 if( Ndb1apd[ia]>=28 ) Napd_28cells[0]+=1.;
1412 if( Ndb2apd[ia]>=24 ) Napd_24cells[1]+=1.;
1413 if( Ndb2apd[ia]>=28 ) Napd_28cells[1]+=1.;
1415 if( Ndb3apd[ia]>=24 ) Napd_24cells[2]+=1.;
1416 if( Ndb3apd[ia]>=28 ) Napd_28cells[2]+=1.;
1419 if( Nmcapd[ia]>=24 ) Napd_24cells[3]+=1.;
1420 if( Nmcapd[ia]>=28 ) Napd_28cells[3]+=1.;
1424 int Napd_24cells_tot=0;
1425 int Napd_28cells_tot=0;
1426 for(
int id=0;
id<4; ++
id ){
1427 Napd_24cells_tot += Napd_24cells[
id];
1428 Napd_28cells_tot += Napd_28cells[
id];
1431 float ratio_outtimehits_deadapds=9999.;
1432 if( Napd_28cells_tot>0. ) ratio_outtimehits_deadapds=ncells_outtime*1.0/(Napd_28cells_tot*1.);
1436 for(
int id=0;
id<ndeaddcms; ++
id )
1452 const time_t timeSec = ts.
timeHigh();
1453 struct tm* timeStruct = localtime(&timeSec);
1479 for(
unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
1517 for(
unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
1530 if (!(fmVertex.isValid()))
continue;
1531 std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
1535 if (vert.size() != 1)
continue;
1547 for(
size_t h = 0;
h < chits->size(); ++
h)
1552 if( iacts.empty() ){
1565 std::vector< art::Ptr<simb::MCFlux> > fluxs = fmf.at(0);
1567 if(fluxs.size() == 1){
1579 if( !(nu.
CCNC()==0) )
continue;
1637 if( mclist->size()>0 && fllist->size()>0 ){
1638 for(
unsigned int i_intx=0; i_intx<mclist->size(); ++i_intx){
1641 TLorentzVector Tslnu(nu.
Nu().
Px(),nu.
Nu().
Py(),
1643 TLorentzVector Tmcnu(mcnu.
Nu().
Px(),mcnu.
Nu().
Py(),
1644 mcnu.
Nu().
Pz(),mcnu.
Nu().
E());
1662 std::vector<const sim::Particle*> AllPi0;
1666 if(genpar->
PdgCode() != 111)
continue;
1667 if( genpar->
Vz()<0. || genpar->
Vz()>1300. )
continue;
1668 if( genpar->
Vx()>200. || genpar->
Vx()<-200. )
continue;
1669 if( genpar->
Vy()>200. || genpar->
Vy()<-200. )
continue;
1673 if( gam1->
PdgCode() != 22 )
continue;
1674 if( gam2->
PdgCode() != 22 )
continue;
1678 if(genpar->
E() > epi0){
1684 if( AllPi0.size()>0 ){
1686 if( AllPi0.size()<20 ) npi0=AllPi0.size();
1688 for(
int i=0;
i<npi0; ++
i ){
1705 Epi0=selectedpi0->
E();
1719 std::vector<art::Ptr<rb::Track>> Ktracks = fmKTrack.at(sliceIdx);
1723 if( !fmKTrack.isValid() )
continue;
1724 if( Ktracks.size() != 1 )
continue;
1725 if( Ktracks[0]->TotalLength()<1300. )
continue;
1726 TVector3
dir = Ktracks[0]->Dir();
1727 double track_theta=dir.Theta();
1728 if(
cos(track_theta)<0.9 )
continue;
1729 TVector3
start = Ktracks[0]->Start();
1730 double track_startz = start.Z();
1731 if( track_startz>30. )
continue;
1734 double longestKTrack = 0;
1735 int longestKTrack_ncells=0;
1736 double longestKTrack_energy=0.;
1737 int longestKTrack_nplanes=0;
1738 float longestKTrack_phi=0.;
1739 float longestKTrack_theta=0.;
1740 float longestKTrack_startx=0.;
1741 float longestKTrack_starty=0.;
1742 float longestKTrack_startz=0.;
1743 float longestKTrack_stopx=0.;
1744 float longestKTrack_stopy=0.;
1745 float longestKTrack_stopz=0.;
1746 int longestKTrack_is3D=0;
1748 if( fmKTrack.isValid() ){
1750 for(
unsigned int n = 0;
n < Ktracks.size(); ++
n){
1754 const double L = Ktracks[
n]->TotalLength();
1755 if(L > longestKTrack){
1757 longestKTrack_ncells=Ktracks[
n]->NCell();
1758 longestKTrack_energy=Ktracks[
n]->TotalGeV();
1759 longestKTrack_nplanes=Ktracks[
n]->ExtentPlane();
1760 TVector3
dir = Ktracks[
n]->Dir();
1761 longestKTrack_phi=dir.Phi();
1762 longestKTrack_theta=dir.Theta();
1763 TVector3
start = Ktracks[
n]->Start();
1764 TVector3 stop = Ktracks[
n]->Stop();
1765 longestKTrack_startx = start.X();
1766 longestKTrack_starty = start.Y();
1767 longestKTrack_startz = start.Z();
1768 longestKTrack_stopx = stop.X();
1769 longestKTrack_stopy = stop.Y();
1770 longestKTrack_stopz = stop.Z();
1771 if( Ktracks[
n]->Is3D() ) longestKTrack_is3D=1;
1779 std::vector<float> tcellX;
1780 std::vector<float> tcellY;
1781 std::vector<float> tcellZ;
1782 std::vector<float> tcellE;
1783 std::vector<int> tcellPlane;
1784 std::vector<int> tcellCell;
1785 std::vector<int> tcellADC;
1786 std::vector<float> tcellPECorr;
1787 std::vector<float> tcellTNS;
1788 std::vector<float> tcellW;
1791 for(
unsigned int icell=0; icell<Ktracks[trkIndex]->NCell(); ++ icell){
1794 if( !rhit.IsCalibrated() )
continue;
1798 tcellX.push_back(rhit.X());
1799 tcellY.push_back(rhit.Y());
1800 tcellZ.push_back(rhit.Z());
1801 tcellE.push_back(rhit.GeV());
1802 tcellADC.push_back(chit->
ADC());
1803 tcellPlane.push_back(chit->
Plane());
1804 tcellCell.push_back(chit->
Cell());
1806 tcellW.push_back(Ktracks[trkIndex]->
W(chit.
get()));
1807 tcellPECorr.push_back(rhit.PECorr());
1811 for(
int ic=0; ic<ncells_trk; ++ic ){
1829 double nCells = slice.
NCell();
1833 double recoEnergy=0.;
1836 for(
unsigned int icell=0; icell<slice.
NCell(); ++ icell){
1840 if( !rhit.IsCalibrated() )
continue;
1841 recoEnergy += rhit.
GeV();
1843 if( rhit.PECorr()>100. && rhit.PECorr()<245. ) nmip++;
1847 if( nCells<20 || nCells>500 )
continue;
1850 if(
fabs(vert[0]->
GetX())>140. )
continue;
1851 if(
fabs(vert[0]->
GetY())>140. )
continue;
1852 if( vert[0]->GetZ()<100. )
continue;
1853 if( vert[0]->GetZ()>700. )
continue;
1863 for(
int id=0;
id<ndeaddcms; ++
id )
1937 lid = elid->Value();
1938 lide = elid->ValueWithE();
1944 std::vector<const slid::ShowerLID* > lids = fmShwLID.at(sliceIdx);
1948 if( nshowers>20 ) nshowers=20;
1950 if( lids.size()>0 ){
1953 std::vector<int> shw_index;
1954 std::vector<float> shwX;
1955 std::vector<float> shwY;
1956 std::vector<float> shwZ;
1957 std::vector<float> shwE;
1958 std::vector<int> shwPlane;
1959 std::vector<int> shwCell;
1960 std::vector<int> shwADC;
1961 std::vector<float> shwPECorr;
1963 for(
int is=0; is<
nshowers; ++is ){
1964 shwLID[is] = lids[is]->Value();
1968 if( shw->Is3D() )
shw3D[is]=1;
1993 shwEnu[is] = shw->CalorimetricEnergy();
1995 for(
unsigned int icell=0; icell<shw->NCell(); ++icell ){
1998 if( !rhit.IsCalibrated() )
continue;
2001 shw_index.push_back(is);
2002 shwX.push_back(rhit.X());
2003 shwY.push_back(rhit.Y());
2004 shwZ.push_back(rhit.Z());
2005 shwE.push_back(rhit.GeV());
2006 shwPlane.push_back(chit->
Plane());
2007 shwCell.push_back(chit->
Cell());
2008 shwADC.push_back(chit->
ADC());
2009 shwPECorr.push_back(rhit.PECorr());
2013 if( shw_ncell>500 ) shw_ncell=500;
2015 for(
int ic=0; ic<shw_ncell; ++ic ){
2042 std::vector<art::Ptr<rb::Prong>> SelectedProng = fmProng3D.at(0);
2046 int nprongs=SelectedProng.size();
2047 if( nprongs>20 ) nprongs=20;
2049 std::vector<double> ProngEnergy;
2050 std::vector<double> ProngEnergyXView;
2051 std::vector<double> ProngEnergyYView;
2052 std::vector<double> ProngPhi;
2054 if( SelectedProng.size()<1 )
continue;
2056 if( SelectedProng.size()>0 ){
2058 int allprong_ncell=0;
2059 std::vector<int> prongIndex;
2060 std::vector<float> pcellX;
2061 std::vector<float> pcellY;
2062 std::vector<float> pcellZ;
2063 std::vector<float> pcellE;
2064 std::vector<int> pcellPlane;
2065 std::vector<int> pcellCell;
2066 std::vector<int> pcellADC;
2067 std::vector<float> pcellPECorr;
2068 std::vector<float> pcellW;
2072 TVector3 dir1 = SelectedProng[
ip]->Dir();
2077 TVector3 start1 = SelectedProng[
ip]->Start();
2084 TVector3 Pstop = start1 + (SelectedProng[
ip]->TotalLength())*dir1;
2100 for(
unsigned int icell=0; icell<SelectedProng[
ip]->NCell(); ++ icell){
2103 double wx=SelectedProng[
ip]->W(chit.
get());
2106 if( !rhit.IsCalibrated() )
continue;
2107 Eprong += rhit.
GeV();
2111 EprongX += rhit.GeV();
2113 pcellW.push_back(wx);
2116 EprongY += rhit.GeV();
2118 pcellW.push_back(wy);
2122 prongIndex.push_back(
ip);
2123 pcellX.push_back(rhit.X());
2124 pcellY.push_back(rhit.Y());
2125 pcellZ.push_back(rhit.Z());
2126 pcellE.push_back(rhit.GeV());
2128 pcellADC.push_back(chit->
ADC());
2129 pcellPlane.push_back(chit->
Plane());
2130 pcellCell.push_back(chit->
Cell());
2131 pcellPECorr.push_back(rhit.PECorr());
2134 ProngEnergy.push_back(Eprong);
2135 ProngEnergyXView.push_back(EprongX);
2136 ProngEnergyYView.push_back(EprongY);
2137 ProngPhi.push_back(dir1.Phi());
2146 if( allprong_ncell>500 ) allprong_ncell=500;
2148 for(
int ic=0; ic<allprong_ncell; ++ic ){
2166 double MaxEnergy1=-999.;
2167 double MaxEnergy2=-999.;
2170 if( MaxEnergy1<ProngEnergy[
i] ){
2171 MaxEnergy1=ProngEnergy[
i];
2186 if(
i==Index1 )
continue;
2187 if( MaxEnergy2<ProngEnergy[
i] ){
2188 MaxEnergy2=ProngEnergy[
i];
2192 double Eden=ProngEnergy[Index1]+ProngEnergy[Index2];
2193 double Enum=ProngEnergy[Index1]-ProngEnergy[Index2];
2194 if( Eden>0.0 ) ratio=Enum/Eden;
2195 dphi=
fabs(ProngPhi[Index1]-ProngPhi[Index2]);
2196 if( dphi>2.*TMath::Pi() ) dphi=dphi-2.*TMath::Pi();
2197 if( dphi>TMath::Pi() ) dphi=2.*TMath::Pi()-dphi;
2207 std::vector<float> p1cellX;
2208 std::vector<float> p1cellY;
2209 std::vector<float> p1cellZ;
2210 std::vector<float> p1cellE;
2211 std::vector<int> p1cellPlane;
2212 std::vector<int> p1cellCell;
2213 std::vector<int> p1cellADC;
2214 std::vector<float> p1cellPE;
2215 std::vector<float> p1cellPECorr;
2216 std::vector<int> p1cellXView;
2217 std::vector<float> p1cellTNS;
2218 std::vector<float> p1cellW;
2220 for(
unsigned int icell=0; icell<SelectedProng[Index1]->NCell(); ++ icell){
2223 double wx=SelectedProng[Index1]->W(chit.
get());
2225 if( !rhit.IsCalibrated() )
continue;
2227 if( rhit.PECorr()>100. && rhit.PECorr()<240. ) ++p1_mipcells;
2230 p1cellW.push_back(wx);
2233 p1cellW.push_back(wy);
2236 p1cellX.push_back(rhit.X());
2237 p1cellY.push_back(rhit.Y());
2238 p1cellZ.push_back(rhit.Z());
2239 p1cellE.push_back(rhit.GeV());
2241 p1cellADC.push_back(chit->
ADC());
2242 p1cellPlane.push_back(chit->
Plane());
2243 p1cellCell.push_back(chit->
Cell());
2244 p1cellPECorr.push_back(rhit.PECorr());
2248 if( p1_ncells>500 ) p1_ncells=500;
2251 for(
int ic=0; ic<p1_ncells; ++ic ){
2268 int ncells_slice=0.;
2269 std::vector<float> slicecellX;
2270 std::vector<float> slicecellY;
2271 std::vector<float> slicecellZ;
2272 std::vector<float> slicecellE;
2273 std::vector<int> slicecellPlane;
2274 std::vector<int> slicecellCell;
2275 std::vector<int> slicecellADC;
2276 std::vector<float> slicecellPE;
2277 std::vector<float> slicecellPECorr;
2278 std::vector<int> slicecellXView;
2279 std::vector<float> slicecellTNS;
2280 std::vector<float> slicecellW;
2281 std::vector<int> slicecellNoise;
2283 for(
unsigned int icell=0; icell<slice.
NCell(); ++ icell){
2288 if( !rhit.IsCalibrated() )
continue;
2292 slicecellW.push_back(slice.
W(chit.
get()));
2294 slicecellX.push_back(rhit.X());
2295 slicecellY.push_back(rhit.Y());
2296 slicecellZ.push_back(rhit.Z());
2297 slicecellE.push_back(rhit.GeV());
2299 slicecellADC.push_back(chit->
ADC());
2300 slicecellPlane.push_back(chit->
Plane());
2301 slicecellCell.push_back(chit->
Cell());
2302 slicecellPECorr.push_back(rhit.PECorr());
2303 slicecellTNS.push_back(chit->
TNS()/1000.);
2306 if( ncells_slice>500 ) ncells_slice=500;
2309 for(
int ic=0; ic<ncells_slice; ++ic ){
2327 std::vector<art::Ptr<rb::Prong>> SelectedProng2D = fmProng2D.at(0);
2331 std::vector<int> Prong2DInXView;
2332 std::vector<double> ProngEnergy2D;
2333 int nprongs2D=SelectedProng2D.size();
2334 if( nprongs2D>20 ) nprongs2D=20;
2338 int allprong2d_ncell=0;
2339 std::vector<int> prong2dIndex;
2340 std::vector<float> p2dcellX;
2341 std::vector<float> p2dcellY;
2342 std::vector<float> p2dcellZ;
2343 std::vector<float> p2dcellE;
2344 std::vector<int> p2dcellPlane;
2345 std::vector<int> p2dcellCell;
2346 std::vector<int> p2dcellADC;
2347 std::vector<float> p2dcellPECorr;
2348 std::vector<float> p2dcellW;
2350 for(
int ip=0;
ip<nprongs2D; ++
ip ){
2352 int prongview=SelectedProng2D[
ip]->View();
2357 TVector3 dir1 = SelectedProng2D[
ip]->Dir();
2361 TVector3 start1 = SelectedProng2D[
ip]->Start();
2366 TVector3 Pstop = start1 + (SelectedProng2D[
ip]->TotalLength())*dir1;
2371 Prong2DInXView.push_back(prongview);
2375 for(
unsigned int icell=0; icell<SelectedProng2D[
ip]->NCell(); ++ icell){
2377 const TVector3
mean = SelectedProng2D[
ip]->MeanXYZ();
2380 if( !rhit.IsCalibrated() )
continue;
2381 Eprong += rhit.
GeV();
2385 prong2dIndex.push_back(
ip);
2386 p2dcellX.push_back(rhit.X());
2387 p2dcellY.push_back(rhit.Y());
2388 p2dcellZ.push_back(rhit.Z());
2389 p2dcellE.push_back(rhit.GeV());
2391 p2dcellADC.push_back(chit->
ADC());
2392 p2dcellPlane.push_back(chit->
Plane());
2393 p2dcellCell.push_back(chit->
Cell());
2394 p2dcellPECorr.push_back(rhit.PECorr());
2398 ProngEnergy2D.push_back(Eprong);
2401 if( allprong2d_ncell>500 ) allprong2d_ncell=500;
2403 for(
int ic=0; ic<allprong2d_ncell; ++ic ){
2419 if( max3D<ProngEnergy[
i] ){
2420 max3D=ProngEnergy[
i];
2426 for(
int i=0;
i<nprongs2D; ++
i ){
2427 if( max2D<ProngEnergy2D[
i] ){
2428 max2D=ProngEnergy2D[
i];
2432 double ebalance2D=1.;
2433 if( nprongs>0 && nprongs2D>0 ){
2436 if( Prong2DInXView[index2D]==0 ){
2437 Eden=ProngEnergy2D[index2D]+ProngEnergyXView[index3D];
2438 Enum=ProngEnergy2D[index2D]-ProngEnergyXView[index3D];
2441 Eden=ProngEnergy2D[index2D]+ProngEnergyYView[index3D];
2442 Enum=ProngEnergy2D[index2D]-ProngEnergyYView[index3D];
2444 if( Eden>0. ) ebalance2D=
fabs(Enum)/Eden;
2461 TH1* hPOT = tfs->
make<TH1F>(
"hTotalPOT",
";; POT", 1, 0, 1);
2463 TH1* hSPILL = tfs->
make<TH1F>(
"hTotalSpill",
";; SPILL", 1, 0, 1);
2470 const double L = 810.;
2471 const double ldm = 2.35e-3;
2472 const double sin22t13 = 0.1;
2473 const double sin22t23 = 1.;
2474 const double ssth23 = 0.5;
double E(const int i=0) const
::xsd::cxx::tree::id< char, ncname > id
int allprong2DcellCell[500]
float allprongcellPECorr[500]
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.
std::string fInstLabel3D
label for prong
float allprong2DcellE[500]
float allprong2DcellY[500]
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
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...
double SimpleOscProb(const simb::MCFlux &flux, const simb::MCNeutrino &nu) const
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
void endSubRun(const art::SubRun &sr)
std::string fKTrackLabel
label for slices
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned short Plane() const
const simb::MCParticle & Nu() const
constexpr std::uint32_t timeHigh() 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)
float allprong2DcellX[500]
A collection of associated CellHits.
float sl_Fouttimehits_noisyapd
std::string fProngLabel
label for kalman tracks
charged current quasi-elastic
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
int allprong2DcellPlane[500]
std::string fVertexLabel
label for generator information
T sqr(T x)
More efficient square function than pow(x,2)
float Fouttimehits_noisyapd
float prongEnergyXView[20]
bool CompareByEnergy(const slid::ShowerLID *a, const slid::ShowerLID *b)
int NumberDaughters() const
Horizontal planes which measure Y.
object containing MC flux information
int Daughter(const int i) const
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
std::string fCellHitLabel
label for vertex
unsigned short Cell() const
int allprong2DcellADC[500]
int InteractionType() const
DEFINE_ART_MODULE(ROCKMRE)
const simb::MCParticle & Lepton() const
const sim::Particle * Primary(const int) const
signed long long int deltaspilltimensec
void push_back(Ptr< U > const &p)
static constexpr double L
T get(std::string const &key) const
std::string flidLabel
label for rvp
double GetY(int dcm=1, int feb=0, int pix=0)
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
std::string fInstLabel2D
instance label for prongs 3D
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.
float trackcellPECorr[500]
int allprongcellPlane[500]
float prong1cellPECorr[500]
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.
std::string fDataSpillLabel
float allprong2DcellPECorr[500]
double Vx(const int i=0) const
T * make(ARGS...args) const
void analyze(const art::Event &evt)
const rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
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 MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) 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 deaddcm_location[14]
void reconfigure(const fhicl::ParameterSet &p)
std::string sl_deaddcm_location[14]
bool IsNoise() const
Is the noise flag set?
float prongEnergyYView[20]
double totgoodpot
normalized by 10^12 POT
std::string frvpLabel
label for cell hits
std::string fMCFluxLabel
instance label for prongs 2D
Event generator information.
double spillpot
POT for spill normalized by 10^12.
double Vy(const int i=0) const
Encapsulate the geometry of one entire detector (near, far, ndos)
double GetX(int ndb=14, int db=1, int feb=0, int pix=0)
int allprongcellCell[500]
int NumberOfPrimaries() const
int sl_deltaspilltimensec
float allprong2DcellZ[500]