21 #include "TMultiLayerPerceptron.h" 53 #include "Utilities/func/MathUtil.h" 54 #include "Utilities/AssociationUtil.h" 55 #include "DAQChannelMap/DAQChannelMapConstants.h" 90 TH1D *fHtPlaneDedxE[4][11][200];
91 TH1D *fHtPlaneDedxG[4][11][200];
92 TH1D *fHtPlaneDedxMu[4][11][200];
93 TH1D *fHtPlaneDedxPi0[4][11][200];
94 TH1D *fHtPlaneDedxHad[4][11][200];
95 TH1D *fHtPlaneDedxP[4][11][200];
96 TH1D *fHtPlaneDedxN[4][11][200];
97 TH1D *fHtPlaneDedxPi[4][11][200];
98 TH1D *fHtPlaneDedxEqe[4][11][200];
99 TH1D *fHtPlaneDedxEres[4][11][200];
100 TH1D *fHtPlaneDedxEdis[4][11][200];
101 TH1D *fHtPlaneDedxEcoh[4][11][200];
102 TH1D *fHtPlaneDedxEsg[4][41][200];
104 unsigned fHtExpPlaneE[4][11];
105 unsigned fHtExpPlaneG[4][11];
106 unsigned fHtExpPlaneMu[4][11];
107 unsigned fHtExpPlanePi0[4][11];
108 unsigned fHtExpPlaneHad[4][11];
109 unsigned fHtExpPlaneP[4][11];
110 unsigned fHtExpPlaneN[4][11];
111 unsigned fHtExpPlanePi[4][11];
112 unsigned fHtExpPlaneEqe[4][11];
113 unsigned fHtExpPlaneEres[4][11];
114 unsigned fHtExpPlaneEdis[4][11];
115 unsigned fHtExpPlaneEcoh[4][11];
116 unsigned fHtExpPlaneEsg[4][41];
161 TH1D *fHtCellTransDedxE[4][11][20];
162 TH1D *fHtCellTransDedxG[4][11][20];
163 TH1D *fHtCellTransDedxMu[4][11][20];
164 TH1D *fHtCellTransDedxPi0[4][11][20];
165 TH1D *fHtCellTransDedxHad[4][11][20];
166 TH1D *fHtCellTransDedxP[4][11][20];
167 TH1D *fHtCellTransDedxN[4][11][20];
168 TH1D *fHtCellTransDedxPi[4][11][20];
169 TH1D *fHtCellTransDedxEqe[4][11][20];
170 TH1D *fHtCellTransDedxEres[4][11][20];
171 TH1D *fHtCellTransDedxEdis[4][11][20];
172 TH1D *fHtCellTransDedxEcoh[4][11][20];
173 TH1D *fHtCellTransDedxEsg[4][41][20];
206 double GetPointDoca(TVector3 vtx, TVector3
start, TVector3 stop);
207 TVector3 GetCellDistToPoint(
unsigned int plane,
unsigned int cell, TVector3 point);
208 double GetCellDistToTrk(
unsigned int plane,
unsigned int cell, TVector3
start, TVector3 stop);
216 TVector3 GetTrkHitPos(
unsigned int plane,
unsigned int cell, TVector3 trkdir, TVector3 trkstart, TVector3 trkstop);
220 unsigned int GetTrkCPlaneCell(
unsigned int plane,
unsigned int i);
224 double GetTrkDistToEdge(
rb::Prong shower);
226 unsigned int ContStartPlane(
rb::Prong shower,
unsigned int startp,
int contp);
232 double ECorrMC(
double gev);
236 double PlaneDedxDataMC(
double gev,
int detid);
237 double CellTransDataMC(
double gev,
int detid);
246 double GetPlaneDedxProb(
rb::Prong shower,
int ireg,
int igev,
unsigned int plane,
double dedx,
const int type);
248 double GetInterPlaneDedxProb(
rb::Prong shower,
unsigned int plane,
double dedx,
const int type, TVector3
pos,
int iReg,
int igev0,
int igev1,
double e0,
double e1);
252 double GetDedxLongLL(
rb::Prong shower,
const int type,
unsigned int startp,
int contp);
254 double GetCellTransDedx(
rb::Prong shower,
unsigned int cell);
255 double GetCellTransDedxProb(
rb::Prong shower,
int ireg,
int igev,
unsigned int cell,
double dedx,
const int type);
256 double GetInterCellTransDedxProb(
rb::Prong shower,
unsigned int cell,
double dedx,
const int type, TVector3
pos,
int iReg,
int igev0,
int igev1,
double e0,
double e1);
261 unsigned int NMIPPlane(
rb::Prong shower);
262 unsigned int NMIPPlane(
rb::Prong shower,
unsigned int startp);
267 TVector3 GetCellNodePos(
unsigned int plane1,
unsigned int cell1,
unsigned int plane2,
unsigned int cell2);
326 std::unique_ptr< std::vector<jmshower::JMShower> >& recoshowercol,
330 std::vector<jmshower::JMShower> ExclShowers(
const art::Event& evt);
334 bool LoadDedxHistograms(
int detid);
335 bool LoadTreeWeights(
int detid);
369 pCellHitPHThresh(50.),
401 produces< std::vector<jmshower::JMShower> >();
402 produces< art::Assns<jmshower::JMShower, rb::Cluster> >();
403 produces< art::Assns<jmshower::JMShower, rb::Prong> >();
511 fnameE +=
"e10kdedxhist_fd_nue_cc.root";
512 fnameG +=
"g10kdedxhist_fd_nue_cc.root";
513 fnameMu +=
"mu10kdedxhist_fd_nue_cc.root";
514 fnamePi0 +=
"pi010kdedxhist_fd_nue_cc.root";
516 fnameP +=
"protondedxhist_fd.root";
517 fnameN +=
"neutrondedxhist_fd.root";
518 fnamePi +=
"pidedxhist_fd.root";
519 fnameEqe +=
"e10kdedxhist_fd_nue_ccqe.root";
520 fnameEres +=
"e10kdedxhist_fd_nue_ccres.root";
521 fnameEdis +=
"e10kdedxhist_fd_nue_ccdis.root";
522 fnameEcoh +=
"e10kdedxhist_fd_nue_cccoh.root";
523 fnameEsg +=
"singleehist_fd.root";
541 fnameE +=
"e10kdedxhist_fd_nue_cc.root";
542 fnameG +=
"g10kdedxhist_fd_nue_cc.root";
543 fnameMu +=
"mu10kdedxhist_fd_nue_cc.root";
544 fnamePi0 +=
"pi010kdedxhist_fd_nue_cc.root";
546 fnameP +=
"protondedxhist_fd.root";
547 fnameN +=
"neutrondedxhist_fd.root";
548 fnamePi +=
"pidedxhist_fd.root";
549 fnameEqe +=
"e10kdedxhist_fd_nue_ccqe.root";
550 fnameEres +=
"e10kdedxhist_fd_nue_ccres.root";
551 fnameEdis +=
"e10kdedxhist_fd_nue_ccdis.root";
552 fnameEcoh +=
"e10kdedxhist_fd_nue_cccoh.root";
553 fnameEsg +=
"singleehist_fd.root";
557 fnameE +=
"e10kdedxhist_fd_nue_cc.root";
558 fnameG +=
"g10kdedxhist_fd_nue_cc.root";
559 fnameMu +=
"mu10kdedxhist_fd_nue_cc.root";
560 fnamePi0 +=
"pi010kdedxhist_fd_nue_cc.root";
562 fnameP +=
"protondedxhist_fd.root";
563 fnameN +=
"neutrondedxhist_fd.root";
564 fnamePi +=
"pidedxhist_fd.root";
565 fnameEqe +=
"e10kdedxhist_fd_nue_ccqe.root";
566 fnameEres +=
"e10kdedxhist_fd_nue_ccres.root";
567 fnameEdis +=
"e10kdedxhist_fd_nue_ccdis.root";
568 fnameEcoh +=
"e10kdedxhist_fd_nue_cccoh.root";
569 fnameEsg +=
"singleehist_fd.root";
574 <<
" Bad detector id = " << detid;
580 fnameE +=
"e10kdedxhist_fd_nue_cc.root";
581 fnameG +=
"g10kdedxhist_fd_nue_cc.root";
582 fnameMu +=
"mu10kdedxhist_fd_nue_cc.root";
583 fnamePi0 +=
"pi010kdedxhist_fd_nue_cc.root";
585 fnameP +=
"protondedxhist_fd.root";
586 fnameN +=
"neutrondedxhist_fd.root";
587 fnamePi +=
"pidedxhist_fd.root";
588 fnameEqe +=
"e10kdedxhist_fd_nue_ccqe.root";
589 fnameEres +=
"e10kdedxhist_fd_nue_ccres.root";
590 fnameEdis +=
"e10kdedxhist_fd_nue_ccdis.root";
591 fnameEcoh +=
"e10kdedxhist_fd_nue_cccoh.root";
592 fnameEsg +=
"singleehist_fd.root";
597 fnameE +=
"e10kdedxhist_fd_nue_cc.root";
598 fnameG +=
"g10kdedxhist_fd_nue_cc.root";
599 fnameMu +=
"mu10kdedxhist_fd_nue_cc.root";
600 fnamePi0 +=
"pi010kdedxhist_fd_nue_cc.root";
602 fnameP +=
"protondedxhist_fd.root";
603 fnameN +=
"neutrondedxhist_fd.root";
604 fnamePi +=
"pidedxhist_fd.root";
605 fnameEqe +=
"e10kdedxhist_fd_nue_ccqe.root";
606 fnameEres +=
"e10kdedxhist_fd_nue_ccres.root";
607 fnameEdis +=
"e10kdedxhist_fd_nue_ccdis.root";
608 fnameEcoh +=
"e10kdedxhist_fd_nue_cccoh.root";
609 fnameEsg +=
"singleehist_fd.root";
614 fnameE +=
"e10kdedxhist_fd_nue_cc.root";
615 fnameG +=
"g10kdedxhist_fd_nue_cc.root";
616 fnameMu +=
"mu10kdedxhist_fd_nue_cc.root";
617 fnamePi0 +=
"pi010kdedxhist_fd_nue_cc.root";
619 fnameP +=
"protondedxhist_fd.root";
620 fnameN +=
"neutrondedxhist_fd.root";
621 fnamePi +=
"pidedxhist_fd.root";
622 fnameEqe +=
"e10kdedxhist_fd_nue_ccqe.root";
623 fnameEres +=
"e10kdedxhist_fd_nue_ccres.root";
624 fnameEdis +=
"e10kdedxhist_fd_nue_ccdis.root";
625 fnameEcoh +=
"e10kdedxhist_fd_nue_cccoh.root";
626 fnameEsg +=
"singleehist_fd.root";
630 <<
" Bad detector id = " << detid;
733 TFile hfileE(fnameE.c_str());
734 TFile hfileG(fnameG.c_str());
735 TFile hfileMu(fnameMu.c_str());
736 TFile hfilePi0(fnamePi0.c_str());
738 TFile hfileP(fnameP.c_str());
739 TFile hfileN(fnameN.c_str());
740 TFile hfilePi(fnamePi.c_str());
741 TFile hfileEqe(fnameEqe.c_str());
742 TFile hfileEres(fnameEres.c_str());
743 TFile hfileEdis(fnameEdis.c_str());
744 TFile hfileEcoh(fnameEcoh.c_str());
745 TFile hfileEsg(fnameEsg.c_str());
774 for(
int r=0;
r<4;
r++){
775 for(
int i=0;
i<11;
i++){
791 for(
int r=0;
r<4;
r++){
792 for(
int i=0;
i<41;
i++){
799 for(
int r=0;
r<4;
r++){
800 for(
int i=0;
i<11;
i++){
801 for(
int l=0;
l<200;
l++){
803 sprintf(cht,
"ht_%d_%d_%d",
r,
i,
l);
906 for(
int r=0;
r<4;
r++){
907 for(
int i=0;
i<41;
i++){
908 for(
int l=0;
l<200;
l++){
910 sprintf(cht,
"ht_%d_%d_%d",
r,
i,
l);
976 for(
int r=0;
r<4;
r++){
977 for(
int i=0;
i<11;
i++){
978 for(
int l=0;
l<20;
l++){
980 sprintf(cht,
"httrans_%d_%d_%d",
r,
i,
l);
997 for(
int r=0;
r<4;
r++){
998 for(
int i=0;
i<41;
i++){
999 for(
int l=0;
l<20;
l++){
1001 sprintf(cht,
"httrans_%d_%d_%d",
r,
i,
l);
1042 fnameWeightsShape +=
"weights_shape_fd_elec.txt";
1043 fnameWeightsShapeElec +=
"weights_shower_fd_elec.txt";
1044 fnameWeightsShapeE +=
"weights_shape_fd_elec_E.txt";
1045 fnameWeightsShapeNoCos +=
"weights_shape_fd_elec_NoCos.txt";
1046 fnameWeightsShapeENoCos +=
"weights_shape_fd_elec_E_NoCos.txt";
1047 fnameWeightsShapeEL +=
"weights_shape_fd_elec_E_EL.txt";
1048 fnameWeightsShapeepi0 +=
"weights_shape_fd_elec_E_epi0.txt";
1049 fnameWeightsShapeepi0EL +=
"weights_shape_fd_elec_E_epi0EL.txt";
1055 fnameWeightsShape +=
"weights_shape_fd_elec.txt";
1056 fnameWeightsShapeElec +=
"weights_shower_fd_elec.txt";
1057 fnameWeightsShapeE +=
"weights_shape_fd_elec_E.txt";
1058 fnameWeightsShapeNoCos +=
"weights_shape_fd_elec_NoCos.txt";
1059 fnameWeightsShapeENoCos +=
"weights_shape_fd_elec_E_NoCos.txt";
1060 fnameWeightsShapeEL +=
"weights_shape_fd_elec_E_EL.txt";
1061 fnameWeightsShapeepi0 +=
"weights_shape_fd_elec_E_epi0.txt";
1062 fnameWeightsShapeepi0EL +=
"weights_shape_fd_elec_E_epi0EL.txt";
1068 fnameWeightsShape +=
"weights_shape_fd_elec.txt";
1069 fnameWeightsShapeElec +=
"weights_shower_fd_elec.txt";
1070 fnameWeightsShapeE +=
"weights_shape_fd_elec_E.txt";
1071 fnameWeightsShapeNoCos +=
"weights_shape_fd_elec_NoCos.txt";
1072 fnameWeightsShapeENoCos +=
"weights_shape_fd_elec_E_NoCos.txt";
1073 fnameWeightsShapeEL +=
"weights_shape_fd_elec_E_EL.txt";
1074 fnameWeightsShapeepi0 +=
"weights_shape_fd_elec_E_epi0.txt";
1075 fnameWeightsShapeepi0EL +=
"weights_shape_fd_elec_E_epi0EL.txt";
1082 <<
" Bad detector id = " << detid;
1117 double vtxcentergev;
1118 double vtxcentergev5cell;
1119 double vtxcentergev10cell;
1120 double vtxcentergev15cell;
1121 double vtxcentergev20cell;
1130 TTree *trainTree_shape =
new TTree(
"trainTree_shape",
"m_trainTree_shape");
1132 trainTree_shape->Branch(
"egLLL", &egLLL,
"egLLL/D");
1133 trainTree_shape->Branch(
"egLLT", &egLLT,
"egLLT/D");
1134 trainTree_shape->Branch(
"emuLLL", &emuLLL,
"emuLLL/D");
1135 trainTree_shape->Branch(
"emuLLT", &emuLLT,
"emuLLT/D");
1136 trainTree_shape->Branch(
"epi0LLL", &epi0LLL,
"epi0LLL/D");
1137 trainTree_shape->Branch(
"epi0LLT", &epi0LLT,
"epi0LLT/D");
1138 trainTree_shape->Branch(
"epLLL", &epLLL,
"epLLL/D");
1139 trainTree_shape->Branch(
"epLLT", &epLLT,
"epLLT/D");
1140 trainTree_shape->Branch(
"enLLL", &enLLL,
"enLLL/D");
1141 trainTree_shape->Branch(
"enLLT", &enLLT,
"enLLT/D");
1142 trainTree_shape->Branch(
"epiLLL", &epiLLL,
"epiLLL/D");
1143 trainTree_shape->Branch(
"epiLLT", &epiLLT,
"epiLLT/D");
1145 trainTree_shape->Branch(
"dedx0", &dedx0,
"dedx0/D");
1146 trainTree_shape->Branch(
"dedx1", &dedx1,
"dedx1/D");
1147 trainTree_shape->Branch(
"dedx2", &dedx2,
"dedx2/D");
1148 trainTree_shape->Branch(
"dedx3", &dedx3,
"dedx3/D");
1149 trainTree_shape->Branch(
"dedx4", &dedx4,
"dedx4/D");
1150 trainTree_shape->Branch(
"dedx5", &dedx5,
"dedx5/D");
1167 trainTree_shape->Branch(
"gap", &gap,
"gap/D");
1168 trainTree_shape->Branch(
"pi0mass", &pi0mass,
"pi0mass/D");
1169 trainTree_shape->Branch(
"vtxgev", &vtxgev,
"vtxgev/D");
1170 trainTree_shape->Branch(
"shE", &shE,
"shE/D");
1171 trainTree_shape->Branch(
"nuE", &nuE,
"nuE/D");
1172 trainTree_shape->Branch(
"costheta", &costheta,
"costheta/D");
1173 trainTree_shape->Branch(
"length", &length,
"length/D");
1174 trainTree_shape->Branch(
"type", &type,
"type/I");
1176 trainTree_shape->Branch(
"vtxcentergev", &vtxcentergev,
"vtxcentergev/D");
1177 trainTree_shape->Branch(
"vtxcentergev5cell", &vtxcentergev5cell,
"vtxcentergev5cell/D");
1178 trainTree_shape->Branch(
"vtxcentergev10cell", &vtxcentergev10cell,
"vtxcentergev10cell/D");
1179 trainTree_shape->Branch(
"vtxcentergev15cell", &vtxcentergev15cell,
"vtxcentergev15cell/D");
1180 trainTree_shape->Branch(
"vtxcentergev20cell", &vtxcentergev20cell,
"vtxcentergev20cell/D");
1183 mlp_shape =
new TMultiLayerPerceptron(
"@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap,@costheta:22:12:6:type", trainTree_shape);
1184 mlp_shape->LoadWeights(fnameWeightsShape.c_str());
1186 mlp_shape_NoCos =
new TMultiLayerPerceptron(
"@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap:22:12:6:type", trainTree_shape);
1190 mlp_shape_elec =
new TMultiLayerPerceptron(
"@egLLL,@egLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@shE,@gap:18:12:6:type", trainTree_shape);
1194 mlp_shape_E =
new TMultiLayerPerceptron(
"@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap,@costheta,@nuE:22:12:6:type", trainTree_shape);
1195 mlp_shape_E->LoadWeights(fnameWeightsShapeE.c_str());
1198 mlp_shape_E_NoCos =
new TMultiLayerPerceptron(
"@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT,@pi0mass,@shE,@vtxgev,@gap,@nuE:22:12:6:type", trainTree_shape);
1202 mlp_shape_EL =
new TMultiLayerPerceptron(
"@egLLL,@egLLT,@emuLLL,@emuLLT,@epi0LLL,@epi0LLT,@epLLL,@epLLT,@enLLL,@enLLT,@epiLLL,@epiLLT:22:12:6:type", trainTree_shape);
1203 mlp_shape_EL->LoadWeights(fnameWeightsShapeEL.c_str());
1206 mlp_shape_epi0 =
new TMultiLayerPerceptron(
"@dedx0,@dedx1,@dedx2,@dedx3,@costheta:25:12:6:type", trainTree_shape);
1210 mlp_shape_epi0EL =
new TMultiLayerPerceptron(
"@dedx0,@dedx1,@dedx2,@dedx3:25:12:6:type", trainTree_shape);
1239 std::unique_ptr< std::vector<jmshower::JMShower> > recoshowercol(
new std::vector<jmshower::JMShower> );
1240 std::unique_ptr< std::vector<jmshower::JMShower> > evtrecoshowercol(
new std::vector<jmshower::JMShower> );
1254 for(
unsigned int i=0;
i<jmshowercol.size();
i++){
1255 evtrecoshowercol->push_back(jmshowercol[
i]);
1257 unsigned int iSlice = jmshowercol[
i].ISlice();
1261 evt.
put(std::move(evtrecoshowercol));
1262 evt.
put(std::move(assns));
1263 evt.
put(std::move(prongassns));
1271 std::unique_ptr< std::vector<jmshower::JMShower> >&recoshowercol,
1285 std::vector<jmshower::JMShower> evtshowercol;
1286 evtshowercol.clear();
1293 if(
pRecoVtx==
true)
if (!(fmv.isValid()))
return evtshowercol;
1294 if (!(fmprong.isValid()))
return evtshowercol;
1295 for(
unsigned int ic = 0; ic < slicecol->size(); ++ic){
1298 if(clust->
IsNoise())
continue;
1299 std::vector<jmshower::JMShower> showercol;
1302 std::vector<art::Ptr<rb::Prong> > prong2d3d = fmprong.at(ic);
1303 std::vector<art::Ptr<rb::Prong> >
track;
1305 for(
unsigned int iPro = 0; iPro < prong2d3d.size(); ++iPro){
1306 if( prong2d3d[iPro]->Is3D())track.push_back(prong2d3d[iPro]);
1313 if (track.size()==0||track.size()>9999){
1319 std::map<geo::OfflineChan, std::vector<int>> hittrkmap;
1320 std::map<geo::OfflineChan, std::vector<int>> hitshmap;
1324 std::vector<const rb::Cluster* > cluster;
1325 std::vector<rb::Prong> vCluster;
1326 std::vector<rb::Prong> vXCluster;
1327 std::vector<rb::Prong> vYCluster;
1328 std::vector<bool> vNantag;
1329 std::vector<double> vtxGeV;
1330 std::vector<double> vtxCenterGeV;
1331 std::vector<double> trkGeV;
1332 std::vector<double> vtxCenterGeV5cell;
1333 std::vector<double> vtxCenterGeV10cell;
1334 std::vector<double> vtxCenterGeV15cell;
1335 std::vector<double> vtxCenterGeV20cell;
1337 double sliceGeV = 0;
1338 double sliceEtot = 0;
1339 double sliceEdge10GeV = 0;
1340 double sliceEdge5GeV = 0;
1341 double sliceEdge2GeV = 0;
1343 std::vector<int> shNCellToEdge;
1344 std::vector<int> shClosetEdge;
1345 std::vector<double> shEdge20GeV;
1346 std::vector<double> shEdge15GeV;
1347 std::vector<double> shEdge10GeV;
1348 std::vector<double> shEdge5GeV;
1349 std::vector<double> shEdge2GeV;
1351 double sliceExclGeV = 0;
1368 cellhitlist.
clear();
1376 shNCellToEdge.clear();
1377 shClosetEdge.clear();
1378 shEdge20GeV.clear();
1379 shEdge15GeV.clear();
1380 shEdge10GeV.clear();
1384 excludeCluster.
Clear();
1388 std::vector<TVector3> evtvtxcol;
1389 std::vector<unsigned int> evtvtxplanecol;
1393 evtvtxplanecol.clear();
1396 Double_t cellXYZ1[3]={0,0,0};
1397 Double_t cellXYZ2[3]={0,0,0};
1406 std::vector<art::Ptr<rb::Vertex> >
v = fmv.at(ic);
1408 for (
unsigned int j=0;
j<v.size(); ++
j) {
1409 TVector3 evtVtx(0,0,0);
1410 unsigned int vtxPlane = 0;
1412 evtVtx.SetX(v[
j]->
GetX());
1413 evtVtx.SetY(v[
j]->
GetY());
1414 evtVtx.SetZ(v[
j]->GetZ());
1419 for(
int trial = 0; trial < 4; ++trial){
1424 double offsetX = 0, offsetY = 0;
1425 if(trial == 1) offsetX = -2;
1426 if(trial == 2) offsetY = -2;
1427 if(trial == 3) offsetX = offsetY = -2;
1441 evtvtxcol.push_back(evtVtx);
1442 evtvtxplanecol.push_back(vtxPlane);
1447 for(
unsigned int i=0;
i<clust->
NXCell();
i++){
1449 hittrkmap[
key].clear();
1450 hitshmap[
key].clear();
1455 for(
unsigned int i=0;
i<track.size();
i++){
1456 for(
unsigned int ii=0;ii<track[
i]->NCell();ii++){
1458 hittrkmap[
key].push_back((
int)
i);
1459 hitshmap[
key].push_back((
int)i);
1466 const double wx = (clust->
NXCell() > 0) ? clust->
W(clust->
XCell(0).
get()) : 0;
1467 const double wy = (clust->
NYCell() > 0) ? clust->
W(clust->
YCell(0).
get()) : 0;
1470 for(
unsigned int i=0;
i<clust->
NXCell();
i++){
1489 double readDist = xyzd1[3];
1491 if(rhit.IsCalibrated()){
1495 double adc = double(chit->
ADC());
1506 double disttoelec = geom->
Plane(chit->
Plane())->Cell(chit->
Cell())->HalfL()+pigtail-wx;
1519 sliceEdge10GeV += gev;
1530 sliceEdge5GeV += gev;
1541 sliceEdge2GeV += gev;
1548 for(
unsigned int i=0;
i<clust->
NYCell();
i++){
1565 double readDist = xyzd1[3];
1567 if(rhit.IsCalibrated()){
1570 double adc = double(chit->
ADC());
1581 double disttoelec = geom->
Plane(chit->
Plane())->Cell(chit->
Cell())->HalfL()+pigtail-wy;
1595 sliceEdge10GeV += gev;
1605 sliceEdge5GeV += gev;
1615 sliceEdge2GeV += gev;
1624 std::vector<TVector3> startcol,endcol,dircol;
1625 std::vector<TVector3> shstartcol,shendcol;
1626 std::vector<unsigned int> startplanecol,stopplanecol;
1627 std::vector<unsigned int> shstartplanecol,shstopplanecol;
1628 std::vector<bool> backtagcol;
1629 std::vector<double> gapcol;
1630 std::vector<double> vtxdocacol;
1634 startplanecol.clear();
1635 stopplanecol.clear();
1641 shstartplanecol.clear();
1642 shstopplanecol.clear();
1648 for(
unsigned int i=0;
i<track.size();
i++){
1658 TVector3 tStart = track[
i]->Start();
1659 bool nantag =
false;
1660 if(
std::isnan(tStart.x())||
std::abs(tStart.x())>10000){tStart.SetX(track[
i]->MinX());nantag =
true;}
1661 if(
std::isnan(tStart.y())||
std::abs(tStart.y())>10000){tStart.SetY(track[
i]->MinY());nantag =
true;}
1662 if(
std::isnan(tStart.z())||
std::abs(tStart.z())>10000){tStart.SetZ(track[
i]->MinZ());nantag =
true;}
1664 if(
std::isnan(tStop.x()||
std::abs(tStop.x())>10000)){tStop.SetX(track[i]->MaxX());nantag =
true;}
1665 if(
std::isnan(tStop.y()||
std::abs(tStop.y())>10000)){tStop.SetY(track[i]->MaxY());nantag =
true;}
1666 if(
std::isnan(tStop.z()||
std::abs(tStop.z())>10000)){tStop.SetZ(track[i]->MaxZ());nantag =
true;}
1667 TVector3 tVect = (tStop-tStart).
Unit();
1675 vCluster.push_back(Clust);
1676 vXCluster.push_back(Clust);
1677 vYCluster.push_back(Clust);
1678 vNantag.push_back(nantag);
1680 vtxGeV.push_back(0);
1681 vtxCenterGeV.push_back(0);
1682 trkGeV.push_back(0);
1683 shNCellToEdge.push_back(9999);
1684 shClosetEdge.push_back(-1);
1685 shEdge20GeV.push_back(0);
1686 shEdge15GeV.push_back(0);
1687 shEdge10GeV.push_back(0);
1688 shEdge5GeV.push_back(0);
1689 shEdge2GeV.push_back(0);
1691 vtxCenterGeV5cell.push_back(0);
1692 vtxCenterGeV10cell.push_back(0);
1693 vtxCenterGeV15cell.push_back(0);
1694 vtxCenterGeV20cell.push_back(0);
1697 TVector3 evtVtx(track[i]->Start()[0],track[i]->Start()[1],track[i]->Start()[2]);
1702 for(
unsigned iv=0;iv<evtvtxcol.size();iv++){
1706 evtVtx = evtvtxcol[iv];
1714 double vtxToStart = (evtVtx-track[
i]->Start()).
Mag();
1716 double vtxIntToStart =
sqrt(vtxToStart*vtxToStart-vtxDoca*vtxDoca);
1717 double vtxIntToStop =
sqrt(vtxToStop*vtxToStop-vtxDoca*vtxDoca);
1718 TVector3 vtxStart = track[
i]->Start();
1723 TVector3 shStart = track[
i]->Start();
1726 unsigned int shStartPlane = track[
i]->MinPlane();
1727 unsigned int shStopPlane = track[
i]->MaxPlane();
1728 double vtxTrkLength = (vtxStop-vtxStart).
Mag();
1729 bool showerVtxBackTag =
false;
1731 if(vtxToStart>vtxToStop) {
1732 showerVtxBackTag =
true;
1734 vtxStop = track[
i]->Start();
1735 double tmp = vtxIntToStart;
1736 vtxIntToStart = vtxIntToStop;
1739 shStop = track[
i]->Start();
1742 if(track[i]->
Dir()[2] < 0.0){
1744 shStartPlane = track[
i]->MaxPlane();
1745 shStopPlane = track[
i]->MinPlane();
1747 if(vtxIntToStop<vtxTrkLength){
1748 vtxIntToStart = -vtxIntToStart;
1751 TVector3 trkVtxVect = (vtxStop-evtVtx).
Unit();
1752 TVector3 shVect = (shStop-shStart).
Unit();
1756 bool vtxnantag =
false;
1757 if(
std::isnan(tStart.x())||
std::abs(tStart.x())>10000){tStart.SetX(track[i]->MinX());vtxnantag =
true;}
1758 if(
std::isnan(tStart.y())||
std::abs(tStart.y())>10000){tStart.SetY(track[i]->MinY());vtxnantag =
true;}
1759 if(
std::isnan(tStart.z())||
std::abs(tStart.z())>10000){tStart.SetZ(track[i]->MinZ());vtxnantag =
true;}
1760 if(
std::isnan(tStop.x()||
std::abs(tStop.x())>10000)){tStop.SetX(track[i]->MaxX());vtxnantag =
true;}
1761 if(
std::isnan(tStop.y()||
std::abs(tStop.y())>10000)){tStop.SetY(track[i]->MaxY());vtxnantag =
true;}
1762 if(
std::isnan(tStop.z()||
std::abs(tStop.z())>10000)){tStop.SetZ(track[i]->MaxZ());vtxnantag =
true;}
1763 if(vtxnantag==
true) shVect = (shStop-shStart).
Unit();
1764 vCluster[
i].SetStart(shStart);
1765 vCluster[
i].SetDir(shVect);
1766 vXCluster[
i].SetStart(shStart);
1767 vXCluster[
i].SetDir(shVect);
1768 vYCluster[
i].SetStart(shStart);
1769 vYCluster[
i].SetDir(shVect);
1774 startcol.push_back(shStart);
1775 endcol.push_back(shStop);
1776 dircol.push_back(shVect);
1777 startplanecol.push_back(shStartPlane);
1778 stopplanecol.push_back(shStopPlane);
1780 backtagcol.push_back(showerVtxBackTag);
1785 double shGap = vtxIntToStart;
1786 if(vtxIntToStart>=0&&vtxIntToStart<vtxDoca)shGap = vtxDoca;
1787 if(vtxIntToStart<0)shGap = vtxDoca;
1790 gapcol.push_back(shGap);
1792 shstartcol.push_back(shStart);
1793 shendcol.push_back(shStop);
1794 shstartplanecol.push_back(shStartPlane);
1795 shstopplanecol.push_back(shStopPlane);
1803 for (
unsigned int i=0;
i<cellhitlist.
size();
i++) {
1805 std::vector<art::Ptr<rb::CellHit> >::iterator itr = cellhitlist.
begin() +
i;
1811 double cellW = 2*geom->
Plane((*itr)->Plane())->Cell((*itr)->Cell())->HalfW();
1814 double thiteweighttotal = 0;
1815 double thiteweighttotalatt = 0;
1816 double thiteweight[100];
1820 for(
unsigned int itrk=0;itrk<track.size();itrk++){
1821 if(vNantag[itrk]==
true)
continue;
1824 for(
unsigned int imap=0;imap<hittrkmap[
key].size();imap++){
1825 if(hittrkmap[
key][imap]==(
int)itrk)trktag=
true;
1828 if(view ==
geo::kX)vXCluster[itrk].Add((*itr));
1829 if(view ==
geo::kY)vYCluster[itrk].Add((*itr));
1838 TVector3 trkdir=dircol[itrk];
1839 TVector3 trkstart=startcol[itrk];
1840 TVector3 trkstop=endcol[itrk];
1843 thiteweight[itrk]=0;
1847 int plane = (
int)(*itr)->Plane()-(
int)startplanecol[itrk];
1848 if(startplanecol[itrk]>stopplanecol[itrk])plane = -
plane;
1850 for(
unsigned int icc=0;icc<track[itrk]->NCell();icc++){
1858 if(plane>=-8&&plane<0){
1862 double w1 = track[itrk]->W((*itr).get());
1865 double readDist = xyzd1[3];
1871 double witrk = track[itrk]->W((*itr).get());
1873 double pigtail = geom->
GetPigtail(geom->
Plane((*itr)->Plane())->Cell((*itr)->Cell())->Id());
1877 cellgev = rhit.
GeV();
1879 double adc = double((*itr)->ADC());
1890 double adc = double((*itr)->ADC());
1902 vtxGeV[itrk]+=cellgev/ecor;
1903 if(celltrkdist<19.5*cellW)vtxCenterGeV[itrk]+=cellgev/ecor;
1904 if(celltrkdist<=5.0*cellW)vtxCenterGeV5cell[itrk]+=cellgev/ecor;
1905 if(celltrkdist<=10.0*cellW)vtxCenterGeV10cell[itrk]+=cellgev/ecor;
1906 if(celltrkdist<=15.0*cellW)vtxCenterGeV15cell[itrk]+=cellgev/ecor;
1907 if(celltrkdist<=20.0*cellW)vtxCenterGeV20cell[itrk]+=cellgev/ecor;
1910 int shnplane =
std::abs((
int)stopplanecol[
i]-(
int)startplanecol[
i])+1;
1912 if(plane>=-1*shnplane&&plane<2*shnplane&&celltrkdist<19.5*cellW){
1918 double w1 = track[itrk]->W((*itr).get());
1921 double readDist = xyzd1[3];
1923 double witrk = track[itrk]->W((*itr).get());
1925 double pigtail = geom->
GetPigtail(geom->
Plane((*itr)->Plane())->Cell((*itr)->Cell())->Id());
1929 cellgev = rhit.
GeV();
1931 double adc = double((*itr)->ADC());
1942 double adc = double((*itr)->ADC());
1963 if(cmd>dx1){cmd=dx1;ce=1;}
1964 if(cmd>dx2){cmd=dx2;ce=2;}
1965 if(cmd>dy1){cmd=dy1;ce=3;}
1966 if(cmd>dy2){cmd=dy2;ce=4;}
1967 if(cmd>dz1){cmd=dz1;ce=5;}
1968 if(cmd>dz2){cmd=dz2;ce=6;}
1969 if(shNCellToEdge[itrk]>cmd){
1970 shNCellToEdge[itrk]=
cmd;
1971 shClosetEdge[itrk]=ce;
1973 if(cmd<2)shEdge2GeV[itrk] += cellgev/ecor;
1974 if(cmd<5)shEdge5GeV[itrk] += cellgev/ecor;
1975 if(cmd<10)shEdge10GeV[itrk] += cellgev/ecor;
1976 if(cmd<15)shEdge15GeV[itrk] += cellgev/ecor;
1977 if(cmd<20)shEdge20GeV[itrk] += cellgev/ecor;
1982 if((*itr)->Plane()<
std::min(startplanecol[itrk],stopplanecol[itrk])||(*itr)->Plane()>
std::max(startplanecol[itrk],stopplanecol[itrk])){
fHitDistMap[
key].push_back(0);
continue;}
1986 double w1 = track[itrk]->W((*itr).get());
1989 double readDist = xyzd1[3];
2004 if(plane>=8&&celltrkdist>19.5*cellW&&trktag==
false)
continue;
2006 if(0<=plane&&plane<8&&celltrkdist>1.5*cellW&&trktag==
false){
2012 double witrk = track[itrk]->W((*itr).get());
2014 double pigtail = geom->
GetPigtail(geom->
Plane((*itr)->Plane())->Cell((*itr)->Cell())->Id());
2018 cellgev = rhit.
GeV();
2020 double adc = double((*itr)->ADC());
2030 double adc = double((*itr)->ADC());
2041 vtxGeV[itrk]+=cellgev/ecor;
2042 if(celltrkdist<19.5*cellW)vtxCenterGeV[itrk]+=cellgev/ecor;
2043 if(celltrkdist<=5.0*cellW)vtxCenterGeV5cell[itrk]+=cellgev/ecor;
2044 if(celltrkdist<=10.0*cellW)vtxCenterGeV10cell[itrk]+=cellgev/ecor;
2045 if(celltrkdist<=15.0*cellW)vtxCenterGeV15cell[itrk]+=cellgev/ecor;
2046 if(celltrkdist<=20.0*cellW)vtxCenterGeV20cell[itrk]+=cellgev/ecor;
2050 if(
pRecluster==
false&&trktag==
false)
continue;
2052 vCluster[itrk].Add((*itr));
2053 bool isinsh =
false;
2054 for(
unsigned int ish=0;ish<hitshmap[
key].size();ish++){
2055 if(hitshmap[
key][ish]==(
int)itrk)isinsh=
true;
2057 if(isinsh==
false) hitshmap[
key].push_back((
int)itrk);
2059 double witrk = track[itrk]->W((*itr).get());
2061 double pigtail = geom->
GetPigtail(geom->
Plane((*itr)->Plane())->Cell((*itr)->Cell())->Id());
2063 double peattitrk = 0;
2066 peattitrk = rhit.
GeV();
2068 double adc = double((*itr)->ADC());
2075 peattitrk = adc/Att;
2079 const double ai = peattitrk/(*itr)->PE();
2080 thiteweight[itrk] = track[itrk]->TotalPE()*
exp(-3.28769
e-01*celltrkdist);
2082 thiteweighttotal += thiteweight[itrk];
2083 thiteweighttotalatt += thiteweight[itrk]/ai;
2086 for(
unsigned int itrk=0;itrk<track.size();itrk++){
2087 if(vNantag[itrk]==
true){
2093 if(thiteweighttotal>0){
2097 double w1 = track[itrk]->W((*itr).get());
2100 double readDist = xyzd1[3];
2105 gev = (*itr)->PE()*(thiteweight[itrk]/thiteweighttotalatt);
2122 double adc = double((*itr)->ADC());
2129 gev = adc*(thiteweight[itrk]/thiteweighttotal)/Att;
2144 for(
unsigned int itrk=0;itrk<track.size();itrk++){
2145 if(vNantag[itrk]==
true)
continue;
2146 for(
unsigned int i=0;
i<track[itrk]->NCell();
i++){
2151 double w1 = track[itrk]->W(trkcell.
get());
2154 double readDist = xyzd1[3];
2158 double witrk = track[itrk]->W(trkcell.
get());
2166 double adc = double(trkcell->
ADC());
2178 double adc = double(trkcell->
ADC());
2195 for(
unsigned int i=0;
i<clust->
NXCell();
i++){
2212 double readDist = xyzd1[3];
2213 if(rhit.IsCalibrated()){
2217 double adc = double(chit->
ADC());
2228 double disttoelec = geom->
Plane(chit->
Plane())->Cell(chit->
Cell())->HalfL()+pigtail-wx;
2231 sliceExclGeV += gev;
2236 for(
unsigned int i=0;
i<clust->
NYCell();
i++){
2253 double readDist = xyzd1[3];
2254 if(rhit.IsCalibrated()){
2258 double adc = double(chit->
ADC());
2268 double disttoelec = geom->
Plane(chit->
Plane())->Cell(chit->
Cell())->HalfL()+pigtail-wy;
2272 sliceExclGeV += gev;
2277 for(
unsigned int i = 0;
i < vCluster.size();
i++){
2278 std::vector<PCluster> showerpcluster;
2279 std::vector<PCluster> trackpcluster;
2280 showerpcluster.clear();
2281 trackpcluster.clear();
2282 if(vNantag[
i]==
true){
2287 if(startplanecol[
i]<=stopplanecol[
i]){
2288 for(
int ip = (
int)startplanecol[i];
ip<=(
int)stopplanecol[i];
ip++){
2291 showerpcluster.push_back(pcluster);
2292 trackpcluster.push_back(pcluster);
2295 for(
int ip = (
int)startplanecol[i];
ip>=(
int)stopplanecol[i];
ip--){
2298 showerpcluster.push_back(pcluster);
2299 trackpcluster.push_back(pcluster);
2302 for (
unsigned int nh=0;nh<vCluster[
i].NCell();nh++) {
2305 if(startplanecol[i]>stopplanecol[i])plane=-
plane;
2306 showerpcluster[
plane].Add((cellhit));
2309 for(
unsigned int imap=0;imap<hittrkmap[
key].size();imap++){
2310 if(hittrkmap[
key][imap]==(
int)
i)trktag=
true;
2312 if(trktag==
true)trackpcluster[
plane].Add((cellhit));
2317 for(
unsigned int i = 0;
i < vCluster.size();
i++){
2318 double startx=-9999,starty=-9999,startz=-9999;
2319 double stopx=-9999,stopy=-9999,stopz=-9999;
2320 bool startxtag=
false,startytag=
false,startztag=
false;
2321 bool stopxtag=
false,stopytag=
false,stopztag=
false;
2322 if(vNantag[
i]==
true){
2323 vtxdocacol.push_back(9999);
2330 if(planecell==9999)
continue;
2331 Double_t cellXYZ[3]={0,0,0};
2333 if(startztag==
false){startz = cellXYZ[2];startztag=
true;}
2334 if(view==
geo::kX&&startxtag==
false){startx = cellXYZ[0];startxtag=
true;}
2335 if(view==
geo::kY&&startytag==
false){starty = cellXYZ[1];startytag=
true;}
2336 if(startxtag==
true&&startytag==
true&&startztag==
true)
break;
2343 if(planecell==9999)
continue;
2344 Double_t cellXYZ[3]={0,0,0};
2346 if(stopztag==
false){stopz = cellXYZ[2];stopztag=
true;}
2347 if(view==
geo::kX&&stopxtag==
false){stopx = cellXYZ[0];stopxtag=
true;}
2348 if(view==
geo::kY&&stopytag==
false){stopy = cellXYZ[1];stopytag=
true;}
2349 if(stopxtag==
true&&stopytag==
true&&stopztag==
true)
break;
2351 double vtxDoca = 9999;
2352 if(!(startx<-9000||starty<-9000||startz<-9000||stopx<-9000||stopy<-9000||stopz<-9000)&&evtvtxcol.size()>0){
2353 TVector3
start(startx,starty,startz);
2354 TVector3 stop(stopx,stopy,stopz);
2357 vtxdocacol.push_back(vtxDoca);
2393 for(
unsigned int i = 0;
i < vCluster.size();
i++){
2394 if(vNantag[
i]==
true)
continue;
2398 if(vCluster[i].NCell()==0){
2400 std::cout<<
"Zero Cell Track! i, NCell, NTrk "<<i<<
" "<<track[
i]->NCell()<<
" "<<track.size()<<
std::endl;
2401 std::cout<<
"Zero Cell Shower! i, NCell, NSh "<<i<<
" "<<vCluster[
i].NCell()<<
" "<<vCluster.size()<<
std::endl;
2410 shower.
SetXNPlane(vXCluster[i].ExtentPlane());
2411 shower.
SetYNPlane(vYCluster[i].ExtentPlane());
2435 for (
unsigned int nh=0;nh<shower.
NCell();nh++) {
2450 shower.
SetGap(gapcol[i]);
2455 if(evtvtxplanecol.size()>0&&contstartp!=0){
2459 contGap=
std::max(gapcol[i],contGap);
2469 if(trkcor<0.0001||trkcor>1000)trkcor = 1;
2480 if(pos[0]>=0&&pos[1]>=0)iReg = 0;
2481 if(pos[0]>=0&&pos[1]<0)iReg = 1;
2482 if(pos[0]<0&&pos[1]>=0)iReg = 2;
2483 if(pos[0]<=0&&pos[1]<=0)iReg = 3;
2487 double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
2488 double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
2490 for(
int ig=0;ig<11;ig++){
2491 epoint[ig] = (elower[ig]+eupper[ig])/2.;
2496 if(gev>epoint[10]){igev0=10;igev1=10;}
2497 else if(gev<epoint[0]){igev0=0;igev1=0;}
2499 for(
int ig=0;ig<10;ig++){
2500 if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
2503 if(igev0>10)igev0=10;
2504 if(igev1>10)igev1=10;
2505 double e0 = gev-epoint[igev0];
2506 double e1 = epoint[igev1] - gev;
2509 double elowerm[41]= {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
2510 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
2511 10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
2512 15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75
2515 double eupperm[41]= { 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
2516 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
2517 10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
2518 15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75, 20.25
2521 for(
int ig=0;ig<41;ig++){
2522 epointm[ig] = (elowerm[ig]+eupperm[ig])/2.;
2527 if(gev>epointm[40]){igevm0=40;igevm1=40;}
2528 else if(gev<epointm[0]){igevm0=0;igevm1=0;}
2530 for(
int ig=0;ig<40;ig++){
2531 if(gev>epointm[ig]&&gev<epointm[ig+1]){igevm0=ig;igevm1=ig+1;}
2534 if(igevm0>40)igevm0=40;
2535 if(igevm1>40)igevm1=40;
2536 double em0 = gev-epointm[igevm0];
2537 double em1 = epointm[igevm1] - gev;
2582 for(
unsigned int tcell = 0; tcell<20; tcell++){
2592 if(r1>1
e-10&&r2>1
e-10){
2594 if(
plane<10)sum10dr+=dr;
2600 std::vector<TVector3> nodetocorecol;
2601 std::vector<geo::OfflineChan> nodexidcol;
2602 std::vector<geo::OfflineChan> nodeyidcol;
2603 std::vector<double> nodetocoreenergycol;
2604 nodetocorecol.clear();
2605 nodetocoreenergycol.clear();
2612 TVector3 sumENodeToCore(0,0,0);
2620 TVector3 iz = shower.
Dir().Unit();
2621 TVector3 ix0(1.0,0.0,0);
2622 TVector3 iy0(0.0,1.0,0);
2623 TVector3 ix = iy0.Cross(iz).Unit();
2624 TVector3 iy = iz.Cross(ix).Unit();
2633 if(view1==view2)
continue;
2636 std::vector<TVector3> planenodecol;
2637 std::vector<double> planenodeenergycol;
2638 planenodecol.clear();
2639 planenodeenergycol.clear();
2645 double nodez = (pcore - shower.
Start()).
Mag();
2646 TVector3 nodetopcore = node - pcore;
2647 nodetopcore.SetZ(nodez);
2648 TVector3 nodetocore = node - core;
2649 nodetocore.SetZ(nodez);
2652 TVector3 xyz = node - core;
2660 double p = shower.
Dir()[0];
2661 double q = shower.
Dir()[1];
2662 double r = shower.
Dir()[2];
2666 TVector3 vt(q*z0-r*y0,r*x0-p*z0,p*y0-q*x0);
2693 double rcn = nodetocore.Perp();
2697 rcn = rcn*nodeenergy;
2699 sumENodeToCore+=nodeenergy*nodetocore;
2700 I11 += nodeenergy*(nodetocore.Y()*nodetocore.Y());
2701 I12 += -nodeenergy*(nodetocore.X()*nodetocore.Y());
2702 I21 += -nodeenergy*(nodetocore.X()*nodetocore.Y());
2703 I22 += nodeenergy*(nodetocore.X()*nodetocore.X());
2705 planenodecol.push_back(node);
2706 planenodeenergycol.push_back(nodeenergy);
2708 nodetocorecol.push_back(nodetocore);
2709 nodetocoreenergycol.push_back(nodeenergy);
2711 nodexidcol.push_back(key1);
2712 nodeyidcol.push_back(key2);
2714 nodexidcol.push_back(key2);
2715 nodeyidcol.push_back(key1);
2721 double Ixx = 0.5*(I11+I22+
sqrt(
pow(I11+I22,2)-4*(I11*I22-I12*I21)));
2722 double Iyy = 0.5*(I11+I22-
sqrt(
pow(I11+I22,2)-4*(I11*I22-I12*I21)));
2724 if(
std::abs(Ixx+Iyy)>1
e-10)asym=(Ixx-Iyy)/(Ixx+Iyy);
2725 double alpha = 3.14159265/2.;
2727 for(
int inode=0;inode<(
int)nodetocorecol.size();inode++){
2728 nodetocorecol[inode].RotateZ(-alpha);
2731 shower.
SetNodes(nodetocorecol, nodetocoreenergycol,nodexidcol,nodeyidcol);
2732 TVector3 ixyz(Ixx,Iyy,0);
2736 shower.
SetNodes(nodetocorecol, nodetocoreenergycol,nodexidcol,nodeyidcol);
2737 TVector3 ixyz(0,0,0);
2799 showercol.push_back(shower);
2802 std::vector<TLorentzVector> showerpGam;
2805 for(
unsigned int i = 0;
i < showercol.size();
i++){
2806 double showereraw = showercol[
i].Energy();
2807 TLorentzVector showerptrk(showercol[
i].
Dir()[0]*showereraw,showercol[
i].
Dir()[1]*showereraw,showercol[
i].
Dir()[2]*showereraw,showereraw);
2808 showerpGam.push_back(showerptrk);
2811 for(
unsigned int i = 0;
i < showercol.size();
i++){
2812 double minpi0mgg = -9999;
2815 for(
unsigned int j=0;
j< showercol.size();
j++){
2817 if(showercol[
j].ISlice()!=showercol[
i].ISlice())
continue;
2818 TLorentzVector p2g = showerpGam[
i] + showerpGam[
j];
2819 double mgg = p2g.M();
2824 shexcle+=showercol[
j].DepositEnergy();
2826 shexcle+=sliceExclGeV;
2827 if(minpi0mgg<0)minpi0mgg=0.;
2828 showercol[
i].SetPi0Mgg(minpi0mgg);
2829 showercol[
i].SetPi0G2ID(minj);
2830 showercol[
i].SetSliceEtot(sliceEtot);
2831 showercol[
i].SetExclEnergy(shexcle);
2832 double elecE = showercol[
i].Energy();
2833 double hadE = showercol[
i].ExclEnergy();
2835 double nuE = (elecE+hadE/(4.00599e-01
2836 -5.08688e+00*TMath::Exp(-1.07930
e+01*hadE-(-1.17127
e+01)*
pow(hadE,2)-1.01452
e+01*
pow(hadE,3)) +4.71995e+00*TMath::Exp(-1.02536
e+01*hadE-(-9.81034
e+00)*
pow(hadE,2)-4.22648
e+00*
pow(hadE,3))));
2838 showercol[
i].SetNueEnergy(nuE);
2844 for(
unsigned int i = 0;
i < showercol.size();
i++){
2845 double Theta = showerpGam[
i].Angle(nuDir);
2846 showercol[
i].SetTheta(Theta);
2849 for(
unsigned int i=0;
i < showercol.size();
i++){
2867 for(
unsigned int i = 0;
i<showercol.size();
i++){
2868 recoshowercol->push_back(showercol[
i]);
2869 evtshowercol.push_back(showercol[i]);
2872 if((
unsigned int)
fprongidx[
i]>=track.size())
continue;
2879 return evtshowercol;
2887 bool tagPreSel =
true;
2889 std::vector<art::Ptr<rb::Prong> >
track;
2892 for(
unsigned int i = 0;
i < trackhandle->size();
i++){
2897 for(
unsigned int itrk=0;itrk<track.size();itrk++){
2898 if(track[itrk]->ExtentPlane()>200)tagPreSel =
false;
2910 double denorm = (start-stop).
Mag();
2911 if(denorm<0.00001){d = (vtx-
start).
Mag();}
else{
2912 d = ((vtx-
start).Cross(vtx-stop)).
Mag()/denorm;
2920 Double_t cellXYZ[3]={0,0,0};
2924 double celldistx = -99;
2925 double celldisty = -99;
2926 double celldistz = -99;
2928 celldistx =
std::abs(cellXYZ[0]-point[0]);
2930 celldistz =
std::abs(cellXYZ[2]-point[2]);
2933 celldisty =
std::abs(cellXYZ[1]-point[1]);
2934 celldistz =
std::abs(cellXYZ[2]-point[2]);
2936 TVector3 celldist(celldistx,celldisty,celldistz);
2942 Double_t cellXYZ[3]={0,0,0};
2947 double x1 = start[0];
2948 double y1 = start[1];
2949 double z1 = start[2];
2950 double p1 = (stop[0] - start[0])/(stop - start).Mag();
2951 double q1 = (stop[1] - start[1])/(stop - start).Mag();
2952 double r1 = (stop[2] - start[2])/(stop - start).Mag();
2954 double x2 = cellXYZ[0];
2955 double y2 = cellXYZ[1];
2956 double z2 = cellXYZ[2];
2983 dDeNorm1[0][0] =
p1;
2984 dDeNorm1[0][1] = q1;
2985 dDeNorm1[1][0] =
p2;
2986 dDeNorm1[1][1] =
q2;
2987 double det1 = dDeNorm1[0][0]*dDeNorm1[1][1] - dDeNorm1[0][1]*dDeNorm1[1][0];
2988 dDeNorm2[0][0] = q1;
2989 dDeNorm2[0][1] = r1;
2990 dDeNorm2[1][0] =
q2;
2991 dDeNorm2[1][1] = r2;
2992 double det2 = dDeNorm2[0][0]*dDeNorm2[1][1] - dDeNorm2[0][1]*dDeNorm2[1][0];
2993 dDeNorm3[0][0] = r1;
2994 dDeNorm3[0][1] =
p1;
2995 dDeNorm3[1][0] = r2;
2996 dDeNorm3[1][1] =
p2;
2997 double det3 = dDeNorm3[0][0]*dDeNorm3[1][1] - dDeNorm3[0][1]*dDeNorm3[1][0];
2998 double det0 = ma*me*mi + mb*mf*mg + mc*md*mh - mc*me*mg - mb*md*mi - ma*mf*mh;
2999 double celldist =
std::abs(det0)/
sqrt(det1*det1+det2*det2+det3*det3);
3007 TVector3 shwdir=shw.
Dir();
3008 TVector3 shwstart=shw.
Start();
3009 TVector3 shwend = shwstart + shw.
TotalLength()*shwdir;
3014 TVector3 shwdir=shw->
Dir();
3015 TVector3 shwstart=shw->
Start();
3016 TVector3 shwend = shwstart + shw->
TotalLength()*shwdir;
3023 Double_t cellXYZ[3]={0,0,0};
3028 TVector3 trkdir=trk->
Dir();
3029 TVector3 trkstart=trk->
Start();
3033 double celldist = 0;
3037 trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3039 trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3045 celldist=
std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3051 celldist=
std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3055 if(trkdtr<0)trkdtr=0;
3057 TVector3 def(-1000,-1000,-1000);
3058 TVector3 trkhitposxcell(trkdtr,trky,celldist);
3059 TVector3 trkhitposycell(trkdtr,trkx,celldist);
3060 if(view ==
geo::kX){
return trkhitposxcell;}
3061 if(view ==
geo::kY){
return trkhitposycell;}
3067 Double_t cellXYZ[3]={0,0,0};
3070 TVector3 trkdir=trk.
Dir();
3071 TVector3 trkstart=trk.
Start();
3075 double celldist = 0;
3078 trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3079 trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3084 celldist=
std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3088 celldist=
std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3092 if(trkdtr<0)trkdtr=0;
3093 TVector3 def(-1000,-1000,-1000);
3094 TVector3 trkhitposxcell(trkdtr,trky,celldist);
3095 TVector3 trkhitposycell(trkdtr,trkx,celldist);
3096 if(view ==
geo::kX){
return trkhitposxcell;}
3097 if(view ==
geo::kY){
return trkhitposycell;}
3102 Double_t cellXYZ[3]={0,0,0};
3107 double celldist = 0;
3110 trkx = trkstart.x()+(trkdir.x()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3111 trky = trkstart.y()+(trkdir.y()/trkdir.z())*(cellXYZ[2]-trkstart.z());
3115 celldist=
std::abs((cellXYZ[0]-trkstart.x())-((cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z())));
3119 celldist=
std::abs((cellXYZ[1]-trkstart.y())-((cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z())));
3123 if(trkdtr<0)trkdtr=0;
3124 TVector3 def(-1000,-1000,-1000);
3125 TVector3 trkhitposxcell(trkdtr,trky,celldist);
3126 TVector3 trkhitposycell(trkdtr,trkx,celldist);
3127 if(view ==
geo::kX){
return trkhitposxcell;}
3128 if(view ==
geo::kY){
return trkhitposycell;}
3134 Double_t cellXYZ[3]={0,0,0};
3137 TVector3 trkdir=trk.
Dir();
3138 TVector3 trkstart=trk.
Start();
3140 double trkplanex=99999,trkplaney=99999, trkplanez=0;
3142 trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3143 trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3145 trkplanex = (trkstart.x()+trkstop.x())/2.0;
3146 trkplaney = (trkstart.y()+trkstop.y())/2.0;
3148 trkplanez = cellXYZ[2];
3149 TVector3 trkplanepos(trkplanex, trkplaney, trkplanez);
3150 return trkplanepos;;
3179 Double_t cellXYZ[3]={0,0,0};
3180 Double_t cellXYZ1[3];
3181 Double_t cellXYZmax[3];
3187 double cellwx = (cellXYZ1[0] - cellXYZ[0] + cellW)/geom->
Plane(plane)->
Ncells();
3188 double cellwy = (cellXYZ1[1] - cellXYZ[1] + cellW)/geom->
Plane(plane)->
Ncells();
3189 TVector3 trkdir=trk->
Dir();
3190 TVector3 trkstart=trk->
Start();
3192 double trkplanex=99999,trkplaney=99999;
3196 trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3197 trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3201 trkplanex = (trkstart.x()+ trkstop.x())/2.0;
3202 trkplaney = (trkstart.y()+ trkstop.y())/2.0;
3215 int trkcellId = -9999;
3219 if(view ==
geo::kX)trkcellId =
int((trkplanex+cellwx/2-cellXYZ[0])/cellwx);
3220 if(view ==
geo::kY)trkcellId =
int((trkplaney+cellwy/2-cellXYZ[1])/cellwy);
3225 return unsigned(trkcellId);
3230 Double_t cellXYZ[3]={0,0,0};
3231 Double_t cellXYZ1[3];
3232 Double_t cellXYZmax[3];
3238 double cellwx = (cellXYZ1[0] - cellXYZ[0] + cellW)/geom->
Plane(plane)->
Ncells();
3239 double cellwy = (cellXYZ1[1] - cellXYZ[1] + cellW)/geom->
Plane(plane)->
Ncells();
3240 TVector3 trkdir=trk.
Dir();
3241 TVector3 trkstart=trk.
Start();
3243 double trkplanex=99999,trkplaney=99999;
3247 trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3248 trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3254 trkplanex = (trkstart.x()+trkstop.x())/2.0;
3255 trkplaney = (trkstart.y()+trkstop.y())/2.0;
3270 int trkcellId = -9999;
3275 if(view ==
geo::kX)trkcellId =
int((trkplanex+cellwx/2-cellXYZ[0])/cellwx);
3276 if(view ==
geo::kY)trkcellId =
int((trkplaney+cellwy/2-cellXYZ[1])/cellwy);
3283 return unsigned(trkcellId);
3298 if(gev<0.00001)
return 9999;
3299 double avgid = sumid/gev;
3303 return (
unsigned int)
cid;
3307 Double_t cellXYZ[3]={0,0,0};
3310 TVector3 trkdir=trk->
Dir();
3311 TVector3 trkstart=trk->
Start();
3313 double trkplanex=99999,trkplaney=99999, trkplanez=0;
3315 trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3316 trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3318 trkplanez = cellXYZ[2];
3320 double xe[4]={9999,9999,9999,9999};
3321 double ye[4]={9999,9999,9999,9999};
3322 double ze[4]={9999,9999,9999,9999};
3339 double xtoedge = 9999.;
3340 double ytoedge = 9999.;
3341 double ztoedge = 9999.;
3343 for(
int ii=0;ii<4;ii++){
3344 if(xtoedge>xe[ii])xtoedge=xe[ii];
3345 if(ytoedge>ye[ii])ytoedge=ye[ii];
3346 if(ztoedge>ze[ii])ztoedge=ze[ii];
3366 TVector3 toedge(xtoedge, ytoedge, ztoedge);
3372 Double_t cellXYZ[3]={0,0,0};
3375 TVector3 trkdir=trk.
Dir();
3376 TVector3 trkstart=trk.
Start();
3378 double trkplanex=99999,trkplaney=99999, trkplanez=0;
3380 trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
3381 trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
3383 trkplanez = cellXYZ[2];
3385 double xe[4]={9999,9999,9999,9999};
3386 double ye[4]={9999,9999,9999,9999};
3387 double ze[4]={9999,9999,9999,9999};
3404 double xtoedge = 9999.;
3405 double ytoedge = 9999.;
3406 double ztoedge = 9999.;
3408 for(
int ii=0;ii<4;ii++){
3409 if(xtoedge>xe[ii])xtoedge=xe[ii];
3410 if(ytoedge>ye[ii])ytoedge=ye[ii];
3411 if(ztoedge>ze[ii])ztoedge=ze[ii];
3434 TVector3 toedge(xtoedge, ytoedge, ztoedge);
3440 int i = shower.
ID();
3441 bool isfiducial =
true;
3445 if(disttoedge.x()>2*planerad&&disttoedge.y()>2*planerad&&disttoedge.z()>2*planerad)
continue;
3452 int i = shower.
ID();
3453 double mindist = 999999.;
3456 double minxyz=999999.;
3457 for(
int ii=0;ii<3;ii++){
3460 if(minxyz<mindist)mindist=minxyz;
3467 int i = shower.
ID();
3468 Double_t cellXYZ[3]={0,0,0};
3471 TVector3 trkdir=shower.
Dir();
3472 TVector3 trkstart=shower.
Start();
3474 double trklength =
std::abs((trkstart - trkstop).
Mag());
3477 if(cellL<cellD||cellL<0.0001){
3487 double trkhitpath=cellD*(
sqrt(trkdir.x()*trkdir.x()+trkdir.y()*trkdir.y()+trkdir.z()*trkdir.z())/
std::abs(trkdir.z()));
3493 int i = shower.
ID();
3495 unsigned int newstartp = startp;
3498 for(
int ic=0;ic<contp;ic++){
3528 for (
unsigned int nh=0;nh<cluster.
NCell();nh++) {
3539 int i = shower.
ID();
3540 for (
unsigned int nh=0;nh<shower.
NCell();nh++) {
3551 for (
unsigned int nh=0;nh<cluster.
NCell();nh++) {
3564 int i = shower.
ID();
3565 for (
unsigned int nh=0;nh<shower.
NCell();nh++) {
3576 double poscor=1.+6.28813e-05*avgx+6.15165e-05*avgy;
3577 if(
pCorAbsCell==
false) poscor = 1.0+1.04866e-04*avgx+6.40668e-05*avgy;
3578 gev = gev/(ecor*poscor);
3586 double ecor = 5.68142e-01;
3593 double Att = 4.08050e+03;
3596 Att = 4.08574e+03+3.49097e+04*TMath::Exp(-2.67297
e-03*d)
3597 -1.62442e+04*TMath::Exp(-3.14363
e-02*d-(-3.93281
e-05)*d*d-1.49793
e-08*d*d*d);
3599 Att = 4.09347e+03+3.50384e+04*TMath::Exp(-2.68025
e-03*d)
3600 -1.74280e+04*TMath::Exp(-3.08694
e-02*d-(-3.86870
e-05)*d*d-1.48737
e-08*d*d*d);
3604 Att = 4.08574e+03+3.49097e+04*TMath::Exp(-2.67297
e-03*d)
3605 -1.62442e+04*TMath::Exp(-3.14363
e-02*d-(-3.93281
e-05)*d*d-1.49793
e-08*d*d*d);
3607 Att = 4.09347e+03+3.50384e+04*TMath::Exp(-2.68025
e-03*d)
3608 -1.74280e+04*TMath::Exp(-3.08694
e-02*d-(-3.86870
e-05)*d*d-1.48737
e-08*d*d*d);
3620 ecor = 1.60823+0.000519042*w-2.47044e-07*w*
w;
3622 ecor = 1.3538+0.00051901*w-1.86311e-07*w*
w;
3626 ecor = 1.60823+0.000519042*w-2.47044e-07*w*
w;
3628 ecor = 1.3538+0.00051901*w-1.86311e-07*w*
w;
3632 if(ecor<0.0001)ecor = 1.0;
3645 if(ecor<0.0001)ecor = 1.0;
3658 if(ecor<0.0001)ecor = 1.0;
3663 int i = shower.
ID();
3666 for (
unsigned int nh=0;nh<shower.
NCell();nh++) {
3675 if(gev>0.001)rad = rad/gev;
3680 int i = shower.
ID();
3683 for (
unsigned int nh=0;nh<cluster.
NCell();nh++) {
3692 if(gev>0.001)rad = rad/gev;
3700 for (
unsigned int nh=0;nh<cluster.
NCell();nh++) {
3703 double cellEnergy = cellhit->
PE();
3705 rad += cellEnergy*
sqrt(celldist[0]*celldist[0]+celldist[1]*celldist[1]+celldist[2]*celldist[2]);
3708 if(gev>0.001)rad = rad/gev;
3716 double xadccellx = 0;
3717 double yadccelly = 0;
3718 double zadccellz = 0;
3723 for(
unsigned int nh=0; nh<cluster.
NXCell();nh++){
3725 Double_t cellXYZ[3]={0,0,0};
3727 geom->
Plane(cellhit->
Plane())->Cell(cellhit->
Cell())->GetCenter(cellXYZ);
3728 double cellADC = double(cellhit->
ADC());
3731 xadccellx += cellADC*cellXYZ[0];
3732 zadccellz += cellADC*cellXYZ[2];
3735 for(
unsigned int nh=0; nh<cluster.
NYCell();nh++){
3737 Double_t cellXYZ[3]={0,0,0};
3739 geom->
Plane(cellhit->
Plane())->Cell(cellhit->
Cell())->GetCenter(cellXYZ);
3740 double cellADC = double(cellhit->
ADC());
3743 yadccelly += cellADC*cellXYZ[1];
3744 zadccellz += cellADC*cellXYZ[2];
3746 if(xadc>0)cx = xadccellx/xadc;
3747 if(yadc>0)cy = yadccelly/yadc;
3748 if(zadc>0)cz = zadccellz/zadc;
3749 TVector3 centroid(cx, cy, cz);
3761 int i = shower.
ID();
3773 double dedx = ep/
path;
3785 if(ireg<0||ireg>3)
return 1;
3786 if(igev>=41||igev<0)
return 1;
3787 if(plane>=200)
return 1;
3789 if(ireg<0||ireg>3)
return 1;
3790 if(igev>=11||igev<0)
return 1;
3791 if(plane>=200)
return 1;
3798 if(n<0.0001)n=0.0001;
3810 if(n<0.0001)n=0.0001;
3821 if(n<0.0001)n=0.0001;
3832 if(n<0.0001)n=0.0001;
3857 if(n<0.0001)n=0.0001;
3869 if(n<0.0001)n=0.0001;
3881 if(n<0.0001)n=0.0001;
3892 if(n<0.0001)n=0.0001;
3904 if(n<0.0001)n=0.0001;
3916 if(n<0.0001)n=0.0001;
3928 if(n<0.0001)n=0.0001;
3940 if(n<0.0001)n=0.0001;
4312 double probPXPY = 0;
4313 bool interPos =
false;
4323 double cos = shower.
Dir()[2];
4331 double probreg[4]={1,1,1,1};
4332 for(
int ireg=0;ireg<4;ireg++){
4333 unsigned int plane0 =
int(
int(plane)/
std::abs(cos));
4334 unsigned int plane1 =
int(
int(plane)/
std::abs(cos)) + 1;
4335 double r0 = double(plane/
std::abs(cos))-double(plane0);
4337 double probe0r0 = 1;
4338 double probe0r1 = 1;
4339 double probe1r0 = 1;
4340 double probe1r1 = 1;
4350 double probe0 = (r1*probe0r0+r0*probe0r1);
4351 double probe1 = (r1*probe1r0+r0*probe1r1);
4352 if(igev0!=igev1)probreg[ireg] = (e1*probe0+e0*probe1)/(e0+e1);
4353 if(igev0==igev1)probreg[ireg] = probe0;
4358 double probX1Y1 = probreg[0];
4359 double probX1Y0 = probreg[1];
4360 double probX0Y1 = probreg[2];
4361 double probX0Y0 = probreg[3];
4363 double probX1PY = ((PY-Y0)/(Y1-Y0))*(probX1Y1-probX1Y0)+probX1Y0;
4364 double probX0PY = ((PY-Y0)/(Y1-Y0))*(probX0Y1-probX0Y0)+probX0Y0;
4365 probPXPY = ((PX-X0)/(X1-X0))*(probX1PY-probX0PY)+probX0PY;
4370 if(interPos==
false||probPXPY<0){
4371 unsigned int plane0 =
int(
int(plane)/
std::abs(cos));
4372 unsigned int plane1 =
int(
int(plane)/
std::abs(cos)) + 1;
4373 double r0 = double(plane/
std::abs(cos))-double(plane0);
4375 double probe0r0 = 1;
4376 double probe0r1 = 1;
4377 double probe1r0 = 1;
4378 double probe1r1 = 1;
4383 double probe0 = (r1*probe0r0+r0*probe0r1);
4384 double probe1 = (r1*probe1r0+r0*probe1r1);
4385 if(igev0!=igev1)probPXPY = (e1*probe0+e0*probe1)/(e0+e1);
4386 if(igev0==igev1)probPXPY = probe0;
4463 int i = shower.
ID();
4466 if(pos[0]>=0&&pos[1]>=0)iReg = 0;
4467 if(pos[0]>=0&&pos[1]<0)iReg = 1;
4468 if(pos[0]<0&&pos[1]>=0)iReg = 2;
4469 if(pos[0]<=0&&pos[1]<=0)iReg = 3;
4472 double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
4473 double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
4475 for(
int ig=0;ig<11;ig++){
4476 epoint[ig] = (elower[ig]+eupper[ig])/2.;
4480 if(gev>epoint[10]){igev0=10;igev1=10;}
4481 else if(gev<epoint[0]){igev0=0;igev1=0;}
4483 for(
int ig=0;ig<10;ig++){
4484 if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
4487 if(igev0>10)igev0=10;
4488 if(igev1>10)igev1=10;
4489 double e0 = gev-epoint[igev0];
4490 double e1 = epoint[igev1] - gev;
4492 double elowerm[41]= {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
4493 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
4494 10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
4495 15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75
4498 double eupperm[41]= { 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75,
4499 5.25, 5.75, 6.25, 6.75, 7.25, 7.75, 8.25, 8.75, 9.25, 9.75,
4500 10.25,10.75,11.25,11.75,12.25,12.75,13.25,13.75,14.25,14.75,
4501 15.25,15.75,16.25,16.75,17.25,17.75,18.25,18.75,19.25,19.75, 20.25
4504 for(
int ig=0;ig<41;ig++){
4505 epointm[ig] = (elowerm[ig]+eupperm[ig])/2.;
4510 if(gev>epointm[40]){igevm0=40;igevm1=40;}
4511 else if(gev<epointm[0]){igevm0=0;igevm1=0;}
4513 for(
int ig=0;ig<40;ig++){
4514 if(gev>epointm[ig]&&gev<epointm[ig+1]){igevm0=ig;igevm1=ig+1;}
4517 if(igevm0>40)igevm0=40;
4518 if(igevm1>40)igevm1=40;
4519 double em0 = gev-epointm[igevm0];
4520 double em1 = epointm[igevm1] - gev;
4523 double sumlgprob = 0;
4526 unsigned int nonzeroplane = 0;
4568 if(prob<0.0001)prob=0.0001;
4569 double lgprob =
log(prob);
4582 if(nonzeroplane>0)sumlgprob = sumlgprob/(double)nonzeroplane;
4584 if(totplane>0)sumlgprob = sumlgprob/(double)totplane;
4588 if(sumlgprob>999)sumlgprob=999;
4589 if(sumlgprob<-999)sumlgprob=-999;
4595 int i = shower.
ID();
4598 if(pos[0]>=0&&pos[1]>=0)iReg = 0;
4599 if(pos[0]>=0&&pos[1]<0)iReg = 1;
4600 if(pos[0]<0&&pos[1]>=0)iReg = 2;
4601 if(pos[0]<=0&&pos[1]<=0)iReg = 3;
4604 double elower[11] = {0.00, 0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75};
4605 double eupper[11] = {0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25};
4607 for(
int ig=0;ig<11;ig++){
4608 epoint[ig] = (elower[ig]+eupper[ig])/2.;
4612 if(gev>epoint[10]){igev0=10;igev1=10;}
4613 else if(gev<epoint[0]){igev0=0;igev1=0;}
4615 for(
int ig=0;ig<10;ig++){
4616 if(gev>epoint[ig]&&gev<epoint[ig+1]){igev0=ig;igev1=ig+1;}
4619 if(igev0>10)igev0=10;
4620 if(igev1>10)igev1=10;
4621 double e0 = gev-epoint[igev0];
4622 double e1 = epoint[igev1] - gev;
4625 double sumlgprob = 0;
4628 unsigned int nonzeroplane = 0;
4664 if(prob<0.0001)prob=0.0001;
4665 double lgprob =
log(prob);
4677 if(nonzeroplane>0)sumlgprob = sumlgprob/(double)nonzeroplane;
4679 if(totplane>0)sumlgprob = sumlgprob/(double)totplane;
4682 if(sumlgprob>999)sumlgprob=999;
4683 if(sumlgprob<-999)sumlgprob=-999;
4738 int i = shower.
ID();
4739 double totpath = 0 ;
4754 if(gev<0.00001)
continue;
4755 double avgPos = sumPos/gev;
4757 double wm = avgPos - (double)cid;
4768 int ic = (
int)cell + 1;
4769 double dph = ((double)ic/dir)+wm;
4770 double dpl = ((double)(ic-1)/
dir)+wm;
4772 double eps1 = dph-(
int)dph;
4773 double eps2 = (
int)dpl+1-dpl;
4780 if((
int)fHitEWeightMap[keye2].
size()>i)e2 = eps2*fHitEWeightMap[keye2][
i];